土体的本构模型和超重力物理模拟1)

陈云敏 马鹏程 唐 耀

(浙江大学超重力研究中心,杭州 310058)

引言

土的变形和强度是土力学中的两个重要基本力学问题.土的变形会导致建构筑物沉降,过大的不均匀沉降则会严重影响建构筑物的正常使用,如图1 所示,墨西哥城某建筑和丽江机场即为典型案例.当土的应力水平达到其剪切强度时,会发生剪切破坏,进而可能导致地基失去承载能力和边坡失稳滑坡.如图2 所示,攀枝花机场边坡滑移造成的经济损失超过3 亿元;更为灾难性的深圳光明新区渣土填埋场大滑坡则造成了下游77 人遇难,33 栋房屋倒塌[1].因此,从提高安全性和保障使用功能的角度,在岩土工程设计时掌握土的变形和强度机制,采用合适的方法对土体沉降和稳定性进行定量分析是至关重要的.

图1 土体变形导致建构筑物沉降影响正常使用Fig.1 The effect on normal use of structure caused by the differential settlement of soils

图2 土体剪切破坏导致边坡失稳流滑Fig.2 The slope sliding caused by the shear failure of soils

数值模拟是分析土体沉降和稳定的重要手段.土的本构模型描述了土体复杂的应力应变规律,为岩土工程数值模拟提供了变形和强度的计算公式[2].土是岩石风化产生的碎散工程材料,因此具有压硬性、摩擦性和剪胀性3 大基本力学特性[3].Terzaghi是土力学的奠基人,他提出了有效应力原理和一维固结理论,并将土的基本特性与经典固体力学(弹性多孔介质理论和刚塑性理论)相结合以解决土的沉降和稳定问题[4],引领了传统土力学的发展.但是,传统土力学将土的强度和变形割裂分析的方法是其固有缺陷,因此,一些学者对将塑性力学应用于岩土本构理论的方法进行了持续探索.剑桥大学的Roscoe和他的同事在1963 年提出了第一个可全面反映饱和正常固结黏土基本力学特性的弹塑性模型——剑桥模型(Cam-clay model)[5];随后Roscoe 和Burland 改进了剑桥模型屈服面,提出了更符合土实际变形规律的修正剑桥模型(modified Cam-clay model)[6].剑桥模型通过临界状态理论首次将土的变形和强度统一起来,能够系统地反映土的基本特性,在国际上得到了广泛的接受和应用[7],也标志着土的本构理论进入了新的发展阶段.此后,用于描述不同类型土在不同条件下的变形和强度特性的弹塑性本构模型层出不穷,成为国内外岩土工程学科中较为活跃的研究领域[2].然而,剑桥模型仅适用于饱和正常固结重塑黏土,因此许多学者致力于剑桥模型的改进和发展以适应更多类型的土、更为普遍的应力状态和更加复杂的应力应变关系[8-12].比如,针对超固结土的先剪缩后剪胀、先硬化后软化等不同于正常固结土的复杂力学特性,姚仰平等提出了应力路径无关的统一硬化(unified hardening,UH)参数,并基于UH 参数和修正剑桥模型建立了超固结土的弹塑性本构模型[2,12-13],合理地描述了超固结土的应力应变关系.

物理模拟试验是分析土体沉降和稳定的另一重要手段.物理模拟根据模型尺寸可分为两类,其一是足尺试验,试验模型具有与实际工程相同的尺寸和材料;其二是缩尺试验,试验模型是由原型结构按一定比例缩小得到的,试验数据根据相似准则转化为原型对应的数据.显然,通过足尺试验可获得具有很高参考价值的数据;然而,实际土工建构筑物的尺寸往往很大,受到试验设备、空间和成本等因素的限制,在很多情况下难以进行足尺试验,缩尺试验通常是更为方便和经济的替代选择.常规缩尺试验是在常重力(重力加速度g)下进行的,其缺陷在于无法复现原型的应力水平,而土的变形和强度特性是与应力水平密切相关的,这使缩尺试验结果与原型实际行为之间产生了差距,降低了试验结果的可靠性[14];而超重力物理模拟则可弥补这一缺陷,将缩尺模型置于土工离心机的实验舱内,离心机高速旋转产生的离心力相当于对模型施加了重力加速度高于g的超重力场,进而增大了模型各处的应力水平.1869 年,法国工程师Edouard Phillips 首先提出了超重力试验的设想[15];20 世纪30 年代,美国和前苏联的学者首次采用离心机开展了土工模型超重力试验[16-17];20世纪70 年代以来,伴随着离心机和机载试验装置建设技术的成熟,超重力物理模拟在国际上逐渐得到了广泛的应用,美国、英国、日本等国家相继建设超大型离心机和适用于不同土工试验的先进机载装置[16].我国超重力物理模拟技术虽然起步较晚,但发展比较迅速,自1983 年在长江科学院建成国内首台自主开发的大型土工离心机以来,已有超过二十家高校和科研单位建成了30 余台离心机[14],相关机载装置和试验技术的研发也实现了快速进步.2019 年,国家发改委正式批复由浙江大学牵头建设国家重大科技基础设施超重力离心模拟与实验装置(Centrifugal Hypergravity and Interdisciplinary Experiment Facility,CHIEF),其建成后将是全世界容量最大、应用范围最广的超重力多学科试验平台[18].

本文首先通过比较土和金属的基本物理特性和力学特性,揭示了传统塑性力学中本构模型无法描述土的应力应变关系的事实.然后总结了传统土力学中基于弹性力学的沉降分析方法和基于刚塑性模型的稳定分析方法中存在的问题,阐明了以剑桥模型、统一硬化模型为代表的现代土力学客观反映了土的压硬性、摩擦性和剪胀性,实现了变形和强度、沉降和稳定的统一,并且指出了未来本构模型的发展趋势为考虑土骨架相变和多场耦合.最后,基于土的本构关系,阐述了土体超重力物理模拟的物理本质且其与数值模拟之间的关系,并通过工程实例得到验证.

1 土的基本特性

自然界中的土是岩石的风化、搬运、沉积和再造而产生的,是一种天然的、多相的、碎散的工程材料.土通常是由土颗粒、水和空气组成的三相混合体,土颗粒以一定的结构形式组成土骨架,骨架内部存在孔隙水和孔隙气.孔隙比e是土的重要状态参数,其定义为土中孔隙体积与土骨架体积之比.孔隙比越小,表明土中孔隙的体积相对土骨架的体积越小,即土越致密;反之,则越松散.对于常规类型的土,一般忽略土颗粒自身的变形,因此土的体积变形主要是由孔隙体积的变化导致的,可以通过孔隙比的变化反映出来.

碎散性是土最重要的基本物理特性,是土的力学性质显著区别于金属、混凝土等固体工程材料的根本原因.材料的力学特性由其物理结构决定,因此碎散性也决定了土的基本力学特性,即压硬性、剪胀性和摩擦性[3],进而决定了土的不同于金属等固体材料的应力应变关系.

(1)压硬性是土的变形模量随着均应力的增大而增大的特性.以等向压缩试验为例,其均应力−体应变(p−εv)关系的示意图见图3(a),随着均应力的增大,土产生体积变形,孔隙体积缩小,变形中绝大部分是不可恢复的塑性变形;随着土变得越来越致密,其体积变形模量明显增大,表现出“越压越硬”的特点,体现出了土的压硬性.然而,金属材料的应力应变关系则明显不同.金属的内部晶格之间的排列非常紧密,在等向压力下不会产生晶格间的滑移,因此产生的塑性体应变非常小.

图3 土的基本特性Fig.3 The essential characteristics of soils

图3 土的基本特性(续)Fig.3 The essential characteristics of soils(continued)

(2)剪胀性是土在剪应力作用下体积发生变化的特性,一般将体积膨胀和收缩统称为剪胀性.纯剪切作用下土和金属的体应变−剪应变(εv–εs)关系的示意图见图3(b),可见,土在纯剪切作用下产生了明显的体积应变和剪应变,而金属仅产生剪应变.土的剪胀性是由剪切过程中土骨架变形、土颗粒相对位置变化导致的,对于比较松散的土(如正常固结土和松砂),土颗粒移动后会填充到原先的孔隙中,使其排列变得更为紧密,在宏观上表现为体积的缩小;对于比较致密的土(如超固结土和密砂),原本排列紧密的土颗粒在发生移动后会产生更大的孔隙,使其排列变得更为松散,在宏观上表现为体积的膨胀.金属材料在剪切作用下只发生晶格的相对滑移,不产生体积变形.

(3)摩擦性是土的剪切强度随剪切面上正应力的增大而增大的特性.如图3(c)所示,Tresca 准则是判别金属材料屈服和破坏的常用判据,即破坏面是最大剪应力面,且剪切强度不会因破坏面上的正应力的变化而变化;然而,土的剪切破坏通常使用Mohr-Coulomb 准则判别,剪切破坏发生在剪应力和正应力之比最大的面上.因此,随着破坏面上正应力的增大,剪切强度也会相应增大,这体现了土的摩擦性.土是颗粒材料,剪切破坏时需要克服破坏面上土颗粒间的摩擦力,这是剪切强度与正应力成正比的根本原因,这同时也表明碎散性的基本物理特性对土体力学特性的决定作用.

土的基本力学特性直接影响应力应变关系,除此之外,还有一些力学特性可以通过改变基本特性的发展规律而间接影响了土的应力应变关系,包括各向异性、应力路径依存性等[3].

上文分析了土和金属在应力应变关系上的明显差异,这意味着,以金属本构关系为基础建立的经典塑性力学无法直接用于描述土的应力应变关系.下面以正常固结土的三轴等p剪切试验为例具体说明,在试验中土的均应力p(分别为c1,c2和c3,c1>c2>c3)保持不变,而剪应力q单调增大.采用经典塑性力学中不考虑硬化的理想弹塑性模型得到的应力应变关系如图4(a)所示,而图4(b)为土真实的应力应变关系.剪切强度qf为剪应力的最大值,土的剪切模量G则定义为q–εs曲线的切线斜率.首先,由于土具有压硬性和摩擦性,当均应力p增大时,土的剪切模量G和剪切强度qf都是增大的,但理想弹塑性的应力应变关系中剪切模量G和剪切强度qf均是定值,与p无关;而且,由于土具有剪胀性,在试验中产生了体积变形,这也是理想弹塑性模型无法反映的.

图4 理想弹塑性应力应变关系与土的应力应变关系的区别Fig.4 The differences between the stress-strain relationships of ideal elastic-plasticity and soils

2 土的本构模型

2.1 传统土力学变形和强度理论

1923 年,Terzaghi 提出了土的有效应力和一维固结理论的概念,标志着土力学成为了一门独立的学科.由Terzaghi 主导建立的,以有效应力原理、基于弹性力学的变形分析理论和基于刚塑性模型的破坏分析理论为主要框架的土力学称为传统土力学或古典土力学[20].

饱和土是土颗粒和孔隙水组成的两相混合物,相应地,土受到的外载荷由土骨架和孔隙水共同承担.Terzaghi 认为,真正影响土的变形和强度的并非是孔隙水承担的那一部分应力(即孔隙水压力),而是总应力与孔隙水压力的差值[4],这个差值被定义为效应力.简而言之,土的变形和强度只取决于有效应力.有效应力是土力学的基础,其重要意义在于提供了一座沟通散粒材料和连续固体材料力学分析的“桥梁”,使连续固体力学中的分析方法能够应用于散粒材料.无特殊说明时,本文中提到的土的应力均为有效应力.

传统土力学的应力和变形问题分析采用弹性力学的分析方法,如地基中应力分析的弹性力学解、弹性地基梁分析理论和固结理论等,均将土体视为理想弹性体.固结是土区别于连续固体材料的变形过程,外加载荷作用在饱和土,全部载荷会首先由孔隙水承担,产生超孔隙水压力,随着孔隙水的排出和超孔隙水压力的消散,土的有效应力才会增长,进而导致土骨架变形.传统土力学中的固结理论包括Terzaghi 一维固结理论[4]和Biot 真三维固结理论[21]等,它们将土视为弹性多孔介质,且假设总应力不变.正常固结土在等向压缩试验中的e–p曲线如图5 所示,随着压力水平提高,土的压缩性逐渐减小,这体现了土的压硬性.e–p曲线又被称为正常压缩线(normal compression line,NCL),固结理论通常根据应力水平和NCL 确定压缩模量以考虑压硬性[22-23].

图5 正常固结土的正常压缩线[5]Fig.5 The normal compression line(NCL)of normally consolidated clay[5]

传统土力学在分析土的破坏问题时采用的基于刚塑性模型的极限平衡理论.土被视为刚塑性体,在载荷水平较低时,土体处于静力平衡状态;当载荷达到临界值,土体从静力平衡状态变为运动的临界状态,即极限平衡状态,载荷的临界值称为极限平衡载荷[20],其物理含义是土的破坏载荷.实际上,传统土力学将土的破坏问题转变为求解刚塑性体极限平衡状态的问题,在分析过程中,土的变形是被忽略的.比如在土体稳定分析中广泛应用的Terzaghi 地基极限承载力理论、边坡稳定分析的滑弧分析法等,均是对土中某一滑动破坏面进行极限平衡分析.考虑到土的摩擦性,极限平衡载荷通常根据Mohr-Coulomb 强度准则确定.

传统土力学用相对简单的理论刻画了土不同于连续固体材料的变形和强度特性,应用比较方便,至今仍是岩土工程设计中主要采用的分析方法.但是,传统土力学也有其固有缺陷.首先,传统土力学将土的变形和破坏作为两个割裂的过程分析,分别采用弹性和刚塑性本构,未建立变形和强度之间的关联,这使其无法描述变形逐渐发展直至破坏的过程;其次,将土简化为理想弹性体或刚塑性体会导致较大的分析误差,传统土力学也不能系统地反映土的压硬性、摩擦性和剪胀性以及它们之间的关系,因此若想准确描述土在复杂工况下的应力应变关系,传统土力学往往是无能为力的.

2.2 剑桥模型(Cam-clay model)

传统土力学中对土体特性的简化描述,使其计算结果常常难以满足定量分析的精度需要,因此很多学者致力于发展土的弹塑性本构模型.Roscoe 等[5]提出的剑桥模型(Cam-clay model)作为第一个能够全面反映饱和正常固结重塑土基本特性的弹塑性模型,具有里程碑式的意义;随后,Roscoe 和Burland[6]将剑桥模型的弹头形屈服面修改为椭球圆形,提出了修正剑桥模型(modified Cam-clay model),对土的应力应变关系具有更好的预测效果,至今仍是在工程界中应用最为广泛的土的弹塑性本构模型.

2.2.1 正常固结土的力学特性

剑桥模型是基于正常固结土(后文中的正常固结土和超固结土均指饱和重塑黏土)的试验结果提出的,所谓正常固结,是指土的当前固结压力与历史上受过的最大固结压力相等.Roscoe 将正常固结土称为“wet clay”,这是相较于超固结土而言的,超固结土在历史上受到过更高的固结压力,因此在相同的应力水平下,超固结土的密实度更高,含水量更低.下面分别讨论正常固结土的压硬性、摩擦性和剪胀性以及它们在修正剑桥模型中的描述方法.

(1)压硬性

如前文所述,NCL 是正常固结土压硬性的体现.如图6 所示,NCL 在e−lnp坐标系下近似为直线,剑桥模型用下列方程描述

式中,N为NCL 的截距,物理含义是p=1 kPa 时土的孔隙比;λ 为NCL 的斜率.相应地,等向压缩试验中产生的体积应变为

式中,e0和p0分别为初始孔隙比和初始均应力.修正剑桥模型即采用式(2)描述土的压硬性.卸载时,土发生回弹,回弹线也近似为直线,因此弹性体应变为

因此可得正常固结土在等向压缩试验中产生的塑性体应变为

式中,cp=(λ−κ)/(1+e0).

图6 正常固结土的正常压缩线和回弹线Fig.6 NCL and swelling line of normally consolidated clay

(2)摩擦性

土的摩擦性体现在抗剪强度qf与均应力p成正比,试验表明,二者大致呈直线关系(见图7(a)).因此,修正剑桥模型用如下表达式描述土的摩擦性

式中,M为极限应力比.在临界状态土力学[27]中,土的破坏状态被称为临界状态(critical state),因此M被称为临界状态应力比.

图7 正常固结土的临界状态线Fig.7 Critical state line(CSL)of normally consolidated clay

土体处于临界状态时应力状态和体积(孔隙比)不发生变化,而剪应变快速发展,表现出塑性流动的状态[24].试验表明,正常固结土的应力状态和体积状态是唯一对应的.将p–q–e空间内描述临界状态中应力状态和体积状态唯一对应关系的曲线称为临界状态线(critical state line,CSL),则图7(a)和图7(b)分别为临界状态线在p–q和e–p坐标系内的投影.在p–q坐标系下,CSL 的投影为式(5)所描述的强度包线;在e–p坐标系下,CSL 的投影为与NCL 平行的曲线.临界状态揭示了土体破坏时的应力状态和变形特性,建立了土体变形和强度之间的联系,是剑桥模型的理论基础.

(3)剪胀性

剪胀性反映的是剪应力q的变化对土体体积的影响.描述土的剪胀性,一般采用如下表达式[25-26]

Roscoe 和Burland 基于能量方程推导了修正模型所采用的剪胀方程[6].土单元体塑性功的一般表达式为

再考察图8 中的两个特殊点,点A处于临界状态,只产生塑性剪应变;点B是等向压缩状态,只产生塑性体应变.因此,点A和点B塑性功的表达式分别为

依据式(8),构造了塑性功的特殊表达式

显然,上式是满足式(8)的.联立式(7)和式(9),可得修正剑桥模型的剪胀方程

试验表明,式(10)所示的剪胀方程与正常固结土的试验数据符合程度较好[28].

图8 修正剑桥模型屈服面示意图Fig.8 Schematic diagram of the yield surface of Cam-clay model

2.2.2 修正剑桥模型的屈服面和硬化规律

材料的屈服准则、流动法则和硬化规律是弹塑性模型的3 个基本要素.加工硬化材料屈服面的一般表达式为

式中,H为硬化参数,决定了后继屈服面的变化规律,即材料的硬化规律.修正剑桥模型采用正交流动法则,塑性应变流动方向与塑性势面正交,假设在应力增量dp和dq作用下产生塑性应变增量为,则它们之间满足如下关系

联立式(10)和式(12),可得微分方程

解这个微分方程并考虑边界条件p=px,q=0,可得

式中,px为屈服面与p轴的交点.式(14)即为基于修正剑桥模型剪胀方程和正交流动法则推导的塑性势函数.修正剑桥模型采用相关联流动法则,屈服函数与塑性势函数相同.修正剑桥模型的屈服面如图8 所示,其形状为椭圆.在屈服函数中,硬化参数是px.px的物理含义是等向压缩时的均应力,根据式(4),正常固结土的px与塑性体应变是唯一对应的,因此塑性体应变也可作为硬化参数.联立式(14)和式(4)并进一步整理可得

剑桥模型具有很高的科学意义和应用价值.剑桥模型理论体系完整,参数数量少且可以通过常规三轴试验确定,能够较好地预测饱和正常固结重塑土的应力应变关系,是在岩土工程领域应用最为广泛的本构模型之一.沈珠江[20]指出,“从1963 年Roscoe 发表著名的剑桥模型开始,土力学的研究脱离了古典理论中线弹性多孔介质和刚塑性介质的框框,开始正视土体的非线性变形特性”,因此剑桥模型可被看作是“现代土力学的开端”.

2.3 统一硬化模型(UH model)

自剑桥模型提出后,很多学者致力于它的发展以扩大适用范围,姚仰平等[2,12-13]提出的统一硬化模型(unified hardening model,UH 模型)为典型范例.UH 模型在保持了剑桥模型特色的同时,将适用范围扩大到超固结土.

2.3.1 超固结土的力学特性

正常固结土和超固结土在三轴剪切试验中典型的应力应变关系如图9 所示.正常固结土在剪切过程中始终是硬化、剪缩的,直至达到塑性流动的临界状态(点3);超固结土虽然最终也会达到临界状态,但在剪切过程中则是先硬化(O→2 段)后软化(2→3段)、先剪缩(O→1 段)后剪胀(1→3 段)的.剑桥模型基于正常固结土的试验结果建立,难以描述超固结土的软化和剪胀特性.下面从土的基本特性出发分析超固结土的力学特性及其在UH 模型中的描述方法.

图9 正常固结土和超固结在三轴剪切试验中的应力应变关系Fig.9 Stress-strain relationships in triaxial tests of normally and over consolidated clay

(1)压硬性

如图10,正常固结土在等向卸载时会产生回弹,成为超固结土,而超固结土在压缩时也会产生塑性变形.再压缩线在e–lnp坐标下近似为一条直线,体现了超固结土在再压缩过程中的压硬性.采用与正常固结土类似的推导方法,可得超固结土在等向压缩中的塑性体积应变为

图10 超固结土再压缩线示意图[2]Fig.10 Schematic diagram of the reloading line of over-consolidated clay[2]

(2)摩擦性

图11 中点A为超固结状态点,其均应力为px,对应的前期固结压力点为点B,其均应力为.超固结土的峰值强度高于正常固结土的强度,因此p–q坐标系内的峰值强度包线(称为Hvorslev 线)位于CSL 上方.Hvorslev 线表明超固结土的峰值强度与均应力成正比,体现了超固结土的摩擦性.如图11,当前的Hvorslev 线则由CSL 和过点B的垂线的交点确定,进而点A对应的峰值强度qf可由Hvorslev 线和过点A的垂线的交点确定.根据几何关系,Mf表达式为

式中,M为临界状态应力比,Mh为Hvorslev 线斜率,R0=.

图11 超固结土峰值强度示意图[13]Fig.11 Schematic diagram of the reloading line of over-consolidated line[13]

(3)剪胀性

UH 模型仍然采用修正剑桥模型的剪胀方程(式(10)).如前文所述,土的剪胀方程决定屈服面的形状,因此UH 模型的屈服面也可用式(14)描述.但是,与正常固结土不同的是,超固结土在硬化阶段即出现先剪缩后剪胀现象,也就是说,在px增大时产生了非单调的变化,因此px和不是唯一对应的,不能继续采用作为硬化参数.

2.3.2 统一硬化模型的屈服面和硬化规律

假设UH 模型的硬化参数为H,对比式(15),则超固结土当前屈服面的表达式为

超固结土的超固结程度随着剪切的进行而不断降低,最终达到临界状态,强度也降低至临界状态应力比M,衰化为正常固结土.为描述剪切过程中超固结土的状态变化,姚仰平等[13]提出了参考屈服面、超固结参数以及当前强度的概念.

参考屈服面被定义为与超固结土当前屈服面塑性应变相等的修正剑桥模型屈服面,故超固结状态点A对应的正常固结状态点B即位于参考屈服面上.如图12,假设A′(p,q)是当前屈服面上一点,将它在参考屈服面上相同应力比的应力点记为,则可得参考屈服面的表达式为

其中,R表征了当前屈服面与当前参考屈服面之间的接近程度,即超固结程度.R越大,超固结程度越低;当R=1 时,超固结土衰化为正常固结土.式(17)中的R0为未开始剪切前,土处于等向压缩状态点A时的超固结参数.此外,将超固结土当前状态对应的峰值强度定义为当前强度My,对照式(17),可得My的表达式为

图12 当前屈服面和参考屈服面[13]Fig.12 Current and reference yield surfaces[13]

为描述超固结土屈服面的演化,姚仰平[2,12]等基于等p剪切和等向加载试验的加载特点和结果提出了能统一反映剪胀和剪缩、硬化和软化的统一硬化参数

对于正常固结土,M=My,统一硬化参数可退化为正常固结土的硬化参数,UH 模型可退化为修正剑桥模型.

下面结合超固结土三轴剪切试验各个阶段的应力应变特点讨论UH 模型屈服面和统一硬化参数的合理性.

(a)O→1 阶段,0 <η≤M,特点是剪缩、硬化,点1 (η=M)为剪缩和剪胀的分界点(又称为特征点).根据式(18),UH 模型椭圆屈服面顶点的应力比为M,当应力比低于M时,根据正交流动法则可判断为剪缩变形,故采用UH 模型屈服面可以描述特征点的应力比和此阶段的剪缩变形;在此应力比范围内,dH>0,故采用统一硬化参数可以描述此阶段的硬化.

(b)1→2 阶段,M<η≤My(此处My为 2 点的当前强度),特点是剪胀、硬化,点 2(η=My)为硬化和软化的分界点.当应力比大于M时,根据正交流动法则,采用UH 模型屈服面可以描述该阶段的剪胀变形;在此应力比范围内,dH>0,故采用统一硬化参数可以描述此阶段的硬化.

(c)2→3 阶段,η ≥My≥M,特点是剪胀、软化,在此阶段当前强度降低,应力比也随之降低但略大于当前强度;最终达到临界状态,应力比和当前强度均降低至M(η=My=M,点3).显然,采用UH 模型屈服面可以描述该阶段的剪胀变形;在此应力比范围内,dH<0,故采用统一硬化参数可以描述此阶段的软化.

根据式(21)～式(23),R,My和H均为p,q和的函数,三者的相互关系见图13.当前强度My与超固结程度R有关,硬化参数H与当前强度My有关,同时超固结程度R又与硬化参数H有关,三者之间存在循环的相关关系,因此在采用UH 模型计算时也需要通过循环迭代才可得到应力应变关系.试验表明,UH 模型很好地刻画了超固结土的应力应变关系[29].

图13 R, My 和H 之间的关系Fig.13 The relationship among R, My and H

2.4 土体本构模型的发展趋势

从将强度和变形分离考虑的传统土力学到以剑桥模型、统一硬化模型为代表的全面反映土体力学特性的现代土力学,土体的本构理论实现了长足的进步.目前的本构模型一般认为土体在受力过程中土骨架的颗粒大小和质量不发生变化,土体变形是由土颗粒排列的变化而引起的;但是,在一些工程问题中,岩土体的土骨架会发生颗粒质量损失和级配变化等相变行为,从而导致土体力学特性的改变.显然,目前的本构模型难以描述这些问题中岩土体的应力应变关系.

土骨架不稳定的岩土体广泛存在于岩土工程的不同领域中,有丰富的工程需求.比如,能源岩土工程领域中,在天然气水合物开采和煤热解的过程中分别存在固相向液、气相和固相向气相转化的过程[30-31];环境岩土工程领域中,随着城市固体废弃物的生化降解,固相会向液、气和溶质相转化[32];交通岩土工程领域中,道砟在外力作用下会发生颗粒破裂,进而改变颗粒级配[33];水利岩土工程领域中,高坝中骨料颗粒也有可能在高围压的作用下破裂而发生级配的改变[34].

土骨架变化会显著改变岩土体的力学特性.首先,固相向其他相的转化可能会使土骨架发生弱化,进而降低土体的强度,并产生较大的变形[30,32];其次,级配变化也可能会显著改变土体的强度和变形特性[35].另外,在相变中通常还伴随着力、水、热、生化等多个场的耦合过程,这使其本构关系进一步复杂化.因此,为满足上述复杂岩土问题的定量分析需求,需要发展能考虑土骨架相变和多场耦合的下一代本构模型.

3 土体的超重力物理模拟

土体的超重力物理模拟是指在超重力场中开展物理模拟试验.如图14 所示,超重力场通过离心机产生,将物理模型置于离心机转臂端部的实验舱内,离心机转臂高速旋转给模型提供超过重力加速度g的离心加速度ng,相当于模型受到了ng的超重力作用.超重力场增大了模型的体积力,进而提高了模型内部的应力水平.离心机旋转时,超重力场的方向是水平的.由于模型自身有一定的尺寸,内部不同位置对应的旋转半径不同,因此受到的离心力大小也会有一定差别;离心机转臂的直径越大,模型尺寸对超重力场的影响越小,即模型中的超重力场越均匀.

图14 浙江大学ZJU-400 离心机Fig.14 ZJU-400 centrifuge

3.1 超重力物理模拟的原理

相较常重力场中的物理模拟,超重力物理模拟主要具有缩尺效应和缩时效应两大优势,能够提升对大尺度、长历时的科学和工程问题开展物理模拟的能力.

(1)缩尺效应

在常重力场中,缩尺模型内部的应力水平低于原型,随着缩尺程度的提高,模型与原型差距会越来越大.如上文所述,由于土具有压硬性、摩擦性和剪胀性,其变形和强度与应力水平密切相关.若无法在模型中复现原型的应力水平,即无法复现原型中土体的应力应变关系,进而使模型试验结果失去科学意义及实用价值.但在超重力场中,若模型的缩尺比例为1/n,对其施加ng的超重力加速度即可将模型内部的应力状态提升至与原型一致的水平,在模型土体中复现原型土体的应力应变关系,这即为超重力的缩尺效应.

图15 具体阐述了土的应力状态对其应力应变关系的影响.原型中某一点土体的状态记为点A,位于NCL 上,其均应力为pA,孔隙比为e,其应力应变关系应该与正常固结土相同,不发生剪胀和软化.但若模型缩尺比例为1/100,模型中对应位置土体的状态点则将平移到点A′,孔隙比仍为e,但此时均应力已经降低至0.01pA,土处于超固结状态,会出现图15(b)所示的剪胀和软化现象,相较原型土体存在显著区别,在这种情况下,缩尺模型试验的结果是不具备可靠性的.如图15(b)所示,若将模型置于100g的超重力场中,相应位置的应力状态则能够提升到与原型相同的点A,从而复现原型中正常固结土的应力应变关系.

图15 缩尺模型与实际工程对比示意图Fig.15 Schematic diagram of the comparison between small-scale model and prototype

(2)缩时效应

超重力场增大了不同相物质之间相对运动的驱动力,缩短相物质迁移的时间,因此,超重力的缩时效应即为通过历时较短的模型试验对原型中长历时的演变过程进行模拟.以污染物迁移击穿防污屏障的问题为例,常规模型试验无法完成数十乃至上百年的击穿过程,但若在ng的超重力场中开展试验,模型尺寸可缩小n倍,同时污染物的迁移速度会提高n倍,相应地污染物迁移时间可缩短n2倍.比如,在100g的重力加速度下,模型1 天的对流扩散过程就相当于原型的27.4 年[36].试验时间的大幅缩短使对污染物击穿防污屏障的全过程开展的物理模拟成为可能.此外,缩时效应能够拓展超重力物理模拟应用的科学领域,使其应用于万年时间尺度的地质演变过程重现、新材料的高效探索和制备等方面.

3.2 物理模拟和数值模拟之间的关系

数值模拟和物理模拟是分析土体沉降和稳定的两个重要手段,二者在工程设计中相互配合,相互验证.

随着计算机技术的发展和算法的革新,数值模拟应用广泛.数值模拟的优势在于成本低、速度快,能够实现在物理模拟中难以实现的工况.然而,鉴于数值模拟的基础是本构模型,其结果的可靠性在很大程度上取决于本构模型的合理性,因此在面对本构模型无法合理描述的土体未知特性时,数值模拟的分析结果可靠性往往不高.物理模拟则可以较好地弥补数值模拟的不足,面对无合适本构模型和存在本构模型演变的问题,可通过其寻找规律,揭示土体的未知特性.比如,目前尚无合适的本构模型描述地震中的砂土液化问题,但Zhou 和Chen 通过超重力振动台试验建立了砂土抗液化强度与剪切波速相关模型[37],并得到了美国地质调查局全球地震液化数据库的验证[38].

因此,在面对无合理的本构模型描述土体未知特性时,可以先通过数值模拟建立对此问题的初步认知,再通过超重力物理模拟研究部分工况的响应并将其结果与数值模拟结果进行对比验证,若能够相互验证,后续可用数值模拟分析其他的工况,而不必开展更多的物理模拟试验;若不能相互验证,则需要调整数值模拟的本构模型,或根据物理模拟结果,提出新的理论,并应用于工程分析中.

4 案例——大直径钢管桩水平受荷特性

为研究大直径钢管桩在水平载荷下的响应,开展了多级水平载荷作用下的超重力物理模拟和数值模拟.大直径钢管桩是浙江省某海域海上风机的单桩基础,如图16 所示,桩直径为5.9 m,钢管桩壁厚为0.09 m,桩长79 m,泥面以上高度为24 m,地基土层为海床软黏土[39].

超重力物理模拟采用ZJU-400 离心机,试验重力加速度100g,模型缩尺1/100.地基土采用马来西亚高岭土,力学性质参数为:比重为2.64,液限为60.7%,塑限为36.5%,渗透系数为2.0×10−8m/s,内摩擦角为23◦,有效应力100 kPa 固结系数为每年40 m2,与海床软黏土接近.模型桩采用6061 铝合金空心管,根据相似比,模型桩入土深度55 cm,加载高度24 cm.采用单自由度加载装置对模型桩进行水平加载,载荷共分为10 级,每级1 kN,最大载荷10 kN.采用激光位移计测量位移,桩身应变片测量弯矩.

图16 海上风机单桩基础Fig.16 Single pile foundation of offshore wind turbine

数值分析采用Abaqus 有限元分析软件,地基土的本构采用UH 模型,模型参数见表1,其中前3 个参数根据文献[40]中马来西亚高岭土的试验结果确定;考虑到实际工程中地基土层的排水条件较差,其体积变形几乎为零,故取泊松比ν 为0.49;N是e–lnp坐标系下正常压缩线的纵截距,可根据地基土的初始孔隙比e0、初始应力状态和初始超固结比确定.模型桩采用线弹性模型,弹性模量取68.9 GPa.数值模拟模型尺寸与原型相同.加载模式与超重力物理模拟相同,载荷共分为10 级,根据相似比,每级100 kN,最大载荷1000 kN.

表1 统一硬化模型参数Table 1 Parameters of UH model

超重力物理模拟和数值模拟结果如图17 所示,由此可见,无论是桩身弯矩还是桩身侧向位移,二者的分析结果均比较相近.超重力物理模拟能够验证数值模拟分析结果,证明数值模型和本构模型的合理性.若后续需分析更多的工况(如不同的桩径、桩长和载荷等),为节省成本和时间,可以主要采用数值模拟.

图17 海上风机单桩基础超重力物理模拟和数值模拟结果Fig.17 Results of hypergravity physical and numerical simulations of the single pile foundation of offshore wind turbine

5 结论与展望

(1)土是碎散的工程材料,具有压硬性、摩擦性和剪胀性,这是土的工程力学特性区别于金属的主要特征.

(2)土体的本构模型发展先后历经了强度和变形分离考虑的传统土力学阶段和以剑桥模型、统一硬化模型为代表的全面反映压缩性、摩擦性和剪胀性的现代土力学阶段,对土体复杂应力应变关系的描述效果有明显的提高.未来本构模型的发展趋势为考虑岩土体在受力过程中土骨架相变和多场耦合,旨在为目前本构模型无法定量分析的能源、交通、环境和水利相关的重大岩土工程问题提供合适的数值分析工具.

(3)超重力物理模拟具有缩尺和缩时效应,不仅可验证土体本构模型的合理性,还可以探索现有本构模型无法描述的土体未知特性,在工程中可与数值模拟联合应用,相互验证.除岩土外,超重力物理模拟还可应用于地质、资源和材料等重要领域.国家重大科技基础设施超重力离心模拟与实验装置作为综合超重力多领域学科实验平台,建成后将全面应用于上述领域.