土的组构各向异性及其本构模拟研究进展

杨仲轩，闫珞洢，余 洋，廖 栋，赵朝发

（浙江大学 土木工程学系/岩土工程计算中心/浙江省地下空间开发工程技术研究中心，浙江 杭州 310058）

0 引 言

土是由土颗粒与孔隙组成的多相介质。天然的土体孔隙中常充满水和空气。按土颗粒特性的不同，土通常被分为砂土、黏土、粉土等。按照孔隙中含水量的填充程度不同，我们把同时含有液态水和空气的土称为非饱和土，把仅充满液态水的土称为饱和土。本文重点介绍饱和土的本构模型及其研究进展。

本构模型[1-3]用来定量描述在外荷载作用下土体的力学响应，是土力学的核心研究内容之一，其理论基础是连续介质力学，通常用应力与应变关系表示。早期土体的本构模型主要借鉴于固体（如金属等）。随着对土体力学特性认识的深入，与金属等固体材料不同，土往往表现出多相性、非均质性、非线性、压敏性、剪胀性、各向异性、应力路径与状态相关性、临界状态等复杂力学行为[1,4-6]。为描述这些特性，众多学者采用不同本构假设，构建了多种土体本构模型。

当前主流的土体本构模型通常基于临界状态理论[5]。该理论的提出被认为是土力学发展的里程碑。本文主要以砂土和黏土这两类典型土体为例，对考虑组构各向异性的本构模型研究进展进行总结归纳，并对面临的挑战进行了展望。

本文首先介绍土的组构各向异性，经典临界状态理论和各向异性临界状态理论，随后分别介绍针对砂土和黏土开发的弹塑性本构模型和亚塑性本构模型，最后分析和总结已有模型的优缺点，对如何进一步发展各向异性土体本构模型提出建议。

1 土的组构各向异性

土的组构（Fabric），也称土的内结构（Internal Structure），是对土体内颗粒排布或孔隙分布的宏观统计，通常用二阶张量表示[5,7-10]。组构张量G通常表示为：

式中：N为微观信息数；uk是第k个微观信息统计量（如接触法向、枝向量、颗粒主轴方向等）；w是权重系数；F被称为组构偏张量，其模F代表组构各向异性程度或大小。

实验[9-12]和离散元数值模拟[8,13-14]研究表明：土体的组构及其演化对土体力学特性有重要影响。在本构模型中描述土体的组构各向异性，需确定其初始值、演化规律及最终值。初始值可以通过实验手段获得，但是其最终值很难测定。其根本原因在于土体试样在变形过大时，会产生不同程度的剪切带，导致应变局部化，影响土体临界状态力学特性。多年来，土体组构各向异性临界状态的唯一性也一直是国际学术界争论的焦点，相关内容见后文详述。

2 临界状态理论

临界状态理论包括经典的由 ROSCOE和SCHOFIELD等[15]建立的临界状态理论（Critical State Theory）以及由LI和DAFALIAS[7]建立的各向异性临界状态理论（Anisotropic Critical State Theory）。

2.1 经典临界状态理论

20世纪50年代末至60年代，以剑桥学派为代表建立的临界状态土力学理论[15-18]，为构建土体弹塑性本构模型提供了理论框架。临界状态土力学理论的建立被视为现代土力学发展的里程碑[4-5,15,19]。临界状态是指土体在剪应力作用下会最终达到理想塑性状态，此时土的体积与应力不变，而剪应变则会持续发展。临界状态的数学表达如下：

式中：p为平均正应力；s为偏应力张量；εv为体积应变；e为偏应变张量。由于临界状态给出了土体破坏的终了状态，因此建立本构模型时，描述土体的状态变量可以逐渐逼近该参考状态，而该状态与应力历史和路径无关。自从剑桥模型[18]将临界状态理论引入弹塑性框架后，越来越多的土体本构模型都将临界状态理论作为模型构建的理论基础[9,20-24]。经典临界状态土力学中给出了土体达到临界状态的充分必要条件：孔隙比e和平均正应力p有唯一关系，即达到临界状态线，且应力比η=q/p（q为偏应力第二不变量）等于临界应力比，仅与应力洛德角（或剪切模式）有关。

值得注意的是，经典临界状态理论完全是基于实验观测，并没有任何理论依据。关于临界状态的一些问题，如临界状态线唯一性，临界状态是否能够达到以及平衡态和流动特性等一直以来争论不断[25-32]，其中临界状态是否唯一则是争论的核心问题。大量的室内试验结果表明在e-p平面内临界状态线可能并不唯一，与加载模式（或应力洛德角）、试样的各向异性程度及试样的应力历史等因素有关。而在室内试验中，达到临界状态时由于端部边界条件影响，应力应变状态并不均匀，并通常伴随着剪切带的发展，使获得的临界状态并不可靠。通常认为试样在三轴拉伸条件下比压缩更难达到临界破坏，因为前者达到临界状态时的应变很大，往往超过了仪器设备的能力[9,33]。另一方面，为了描述组构各向异性和剪切模式对土体力学特性的影响，通常引入多条临界状态线作为参考状态，构建相应的本构模型[9,34-38]。需要指出的是，这种处理方式仅仅是本构模拟技巧，并不表明临界状态线不唯一。

2.2 各向异性临界状态理论

经典临界状态理论关于达到临界状态的描述中并没有材料组构各向异性的描述，而基于微观力学的实验和离散元数值模拟表明在临界状态时，土体内结构呈现显著的各向异性，同时土体的力学响应在达到临界状态前也表现为高度各向异性。因此，一种很自然的想法就是在临界状态的描述中引入对组构各向异性的描述。LI和DAFALIAS[7]从热力学基本理论出发，对临界状态问题进行了深入研究。根据吉布斯热力学平衡条件证明了临界状态的唯一性，并且解释了为什么经典临界状态理论中没有包含组构各向异性但并不影响临界状态和临界状态线的原因。尽管如此，考虑到土在临界状态时内结构高度各向异性，组构各向异性是否达到其临界状态值也可以用作判断是否到达临界状态的一个条件，与式（2）同时存在。特别地，如果用一个二阶偏张量来描述土体内结构，在剪切荷载作用下，该张量的模和方向也相应地发生演化。当到达临界状态时，该张量的方向与加载方向一致，而模则达到与应力洛德角相关的临界值。上述条件与经典临界状态理论中关于临界应力比和临界状态线的条件共同构成了临界状态的充分必要条件，作为经典临界状态理论的修正，称为各向异性临界状态理论[7]。

为验证各向异性临界状态理论，YANG和WU[8]开展了离散元数值模拟实验。他们将砂土颗粒近似为重叠小球，生成了不同初始各向异性和孔隙比的试样，获得了在临界状态时砂土的组构各向异性，与各向异性临界状态理论的预测结果吻合。不同应力路径、颗粒形状的离散元数值模拟[13-14,32,39-42]结果也证实了这一结论。ZHAO等[10]和DRUCKREY等[11]通过CT扫描及图像处理技术获得剪切带内砂土颗粒的位移，通过细微观力学统计组构张量、应变与应力的变化，获得了不同密实度砂土组构的临界状态，实验结果验证了各向异性临界状态理论的正确性。

3 各向异性砂土本构模拟

在各向异性临界状态理论框架下，各国学者建立了能够考虑砂土内结构演化的本构模型。其中主要包括弹塑性本构模型和亚塑性本构模型。考虑到砂土的非共轴变形源于其组构各向异性，本文也特别介绍在砂土非共轴模拟方面的研究进展。

3.1 弹塑性本构模型

作为示范性模型，LI和DAFALIAS[7]在土体剪胀和塑性模量方程中引入了组构各向异性变量（组构偏张量与加载方向的第一联合不变量），模拟了三轴空间内主应力方向对初始各向异性砂土的影响。采用相同思路，DAFALIAS及其合作者[43-44]在SANISAND模型框架下，建立了剪胀函数和塑性模量与组构各向异性变量的关系，实现了砂土在不同密度、加载方向和排水条件下的统一模拟。同时，该模型能够模拟主应力轴连续旋转条件下砂土在不同排水条件下的响应。同样，YANG等[45]提出了符合离散元模拟结果的组构演化机制，建立了能考虑不同固结历史对砂土和黏土力学特性影响的统一本构模型。ZHAO及其合作者[46-49]通过在屈服函数和塑性势函数中直接引入组构各向异性变量，用于反映组构各向异性对屈服特性、塑性流动及硬化规律的影响，分别构建了模拟砂土在非三轴空间内的单调加载特性[46,49]、单调和循环荷载[48]和主应力轴旋转路径下[47]的本构模型。

除上述基于经典弹塑性理论的本构模型外，一些学者[50-59]在各自已有模型框架下，采用符合各向异性临界状态理论的组构演化机制，模拟砂土的各向异性力学特性。值得指出的是，ZHAO和KRUYT[59]发展了基于微观力学的满足各向异性临界状态理论的本构模型。与其他模型相比，该模型将组构张量引入在微观量的分布统计函数中，而不是在屈服函数中。

3.2 亚塑性本构模型

亚塑性（Hypoplasticity）模型是一种模拟土体力学特性的本构理论。与弹塑性理论相比，该理论不需把变形分解为弹性和塑性分量，不需引入屈服面、塑性势和硬化准则，仅用一个张量方程描述土体的力学特性。和弹塑性模型相比，亚塑性模型参数相对较少，形式简单，且不需借助相对复杂的塑性力学概念，因此在很多方面具有一定的优越性。

KOLYMBAS[60]创建了亚塑性理论的基本框架，并在该框架下给出了一个示范性的非线性张量函数。其中应力率分解成线性项与非线性项，分别反映土体在加载过程中产生的可恢复和不可恢复变形。WU和BAUER[61]结合砂土的变形特性，提出了第一个较为实用的亚塑性模型，简称为 WU-BAUER亚塑性模型。该模型能较好地模拟土体的剪胀剪缩体积响应以及非线性应力应变关系，且其形式简单，仅需4个模型参数。然而，WU-BAUER亚塑性模型并非基于临界状态理论框架，无法采用同一套参数同时模拟密砂和松砂。另外，由于该模型的剪胀角为常数，不满足土体在临界状态时体应变为零的要求。为解决上述缺陷，VON WOLFFERSDORFF[62]通过在亚塑性模型的非线性项引入孔隙率及临界孔隙率的函数，满足临界状态要求，实现了用同一套参数模拟密砂的剪胀与应变软化，松砂的剪缩与应变硬化等力学行为。为更好地模拟土体在临界状态时的应力状态，VON WOLFFERSDORFF[62]将SMP破坏准则对应的应力状态作为预先定义的临界应力面引入到亚塑性模型中。该模型包含8个参数，且其物理意义较为明确，定量模拟能力获得了较大提升。同样地，GUDEHUS[63]采用三角函数模拟土体的应力状态，提出了8参数的亚塑性本构模型。BAUER[64]通过与实验结果对比，标定了该模型参数，其结果表明 GUDEHUS[63]亚塑性模型拥有 VONWOLFFERSDORFF亚塑性模型类似的模拟能力。

为了模拟各向异性砂土的力学响应，WU[65]在亚塑性非线性项引入各向异性矩阵，反映砂土在不同加载方向上力学行为的差异。然而，该模型并未考虑组构各向异性的演化，且不满足临界状态唯一性的要求。为了解决上述缺陷，合理地模拟组构各向异性及其演化，YANG等[66]将各向异性状态参数引入到VON WOLFFERSDORFF亚塑性模型[62]的特征孔隙率及临界应力面中。该模型满足各向异性临界状态理论[47]的要求。随着加载方向与土体沉积方向趋于一致，各向异性状态参数增加，以模拟由于组构各向异性引起的强度和剪胀的增加。为了模拟在循环荷载作用下各向异性土体的力学特性，LIAO和YANG[52]在YANG等[66]模型中引入了随应变累积而逐渐增大的变量，反映土体在加载过程中的模量衰减。该模型能成功地模拟不排水循环荷载作用下土体应变幅值不断上升、有效应力逐渐减小并最终液化的现象。为模拟土体的非共轴变形，LIAO和YANG[53]将组构张量引入到YANG等[66]的各向异性亚塑性模型的流动方向中，因此流动方向不与应力张量共轴。各向异性临界状态理论[47]在亚塑性框架内的成功应用表明了该理论在各向异性土体本构模拟中的有效性和广泛的适用性。

3.3 非共轴变形模拟

在实际工程问题中，土体的变形往往会伴随着主应力轴的旋转，例如波浪荷载对海床土体的作用，土体遭受多向地震荷载作用，交通荷载作用下路基的变形等。已有研究表明[47,53]，在比例加载条件下，土体非共轴变形的塑性应变较小，而非比例加载产生的非共轴变形则非常明显，不能忽略。在场地条件下，土体通常遭受非比例加载应力路径。因此，土体在非比例加载条件下的非共轴变形是当前研究的热点和难点。

非共轴变形源自土体的组构各向异性。LI和DAFALIAS[34]建立的各向异性临界理论为构建非共轴变形模型提供了较好的理论框架。在临界状态理论框架内，众多学者在非共轴流动和变形模拟方面开展了一系列的研究。

HASHIGUCHI等[67]在次加载面模型中考虑了切向应力增量的影响，通过引入各向异性屈服面来考虑砂土的各向异性并模拟非共轴流动。为建立可以反映主应力轴旋转对砂土变形行为影响的本构模型，GUTIERREZ等[68]利用零弹性区的概念[69]、边界面塑性和亚塑性理论[70]（此处亚塑性与KOLYMBAS[60]建立的亚塑性模型在名称上相同，但其模型框架则完全不同），模拟了应力主轴旋转引起的应力方向和塑性应变增量方向之间的非共轴对土体剪胀的影响，并在剪胀表达式中引入了非共轴系数。虽然该模型建立在二维空间，但在后续工作中通过采用应力和塑性应变增量不变量表示耗散假设，可推广到广义应力空间[71]。LI和DAFALIAS[72]通过切向加载[73]，在LI的边界面模型[74]基础上提出了考虑非比例加载的边界面模型，以反映主应力方向旋转对土体剪切行为的影响。该切向加载机制的剪胀方程和塑性模量都与组构各向异性状态变量相关。值得注意的是，该切向加载机制只有在非比例加载路径下才被激活。LASHKARI和LATIFI[75-76]采用上述非共轴系数引入用于计算非共轴塑性应变的剪胀公式和塑性模量，建立了可以反映纯主应力轴旋转条件下的各向异性的砂土非共轴模型。YANG和 YU[77]基于DAFALIAS和 MANZARI[78-79]的运动硬化及状态相关剪胀的双边界面模型，将塑性应变增量分为由径向应力增量诱导（固定主应力轴，改变主应力大小）和由纯主应力轴旋转引起（应力主值不变，主应力轴方向持续旋转）。每个加载机制都有相应的剪胀方程和塑性模量表达。

4 各向异性黏土本构模拟

黏土颗粒呈片状结构，与砂土颗粒明显不同，因此黏土表现出与砂土不同的物理力学特性。在弹塑性和亚塑性理论框架下，各向异性黏土本构模型也层出不穷。

4.1 弹塑性本构模型

由于沉积过程、固结历史等影响，黏土具有复杂的力学特性，通常表现为不同程度的超固结性和各向异性。为模拟饱和正常固结黏土的力学特性，SCHOFIELD和WROTH[18]建立了首个基于经典临界状态理论的剑桥模型（Cam-Clay Model）。随后，ROSCOE和BURLAND[17]基于更合理的能量耗散假设，建立了修正剑桥模型（Modified Cam-Clay Model）。该模型能更合理地模拟正常固结黏土的压硬性和剪缩特性。与正常固结土相比，超固结土具有较低的孔隙率和较高的强度，随着超固结度的增加，表现出更明显的剪胀和软化特性。在弹塑性理论框架下，通过引入边界面概念[70,80-81]或多重屈服面[82-83]来模拟超固结黏土在初始屈服面内的塑性变形。为描述强超固结黏土的剪胀特性，CHEN和YANG[84-85]在边界面框架下将临界状态线与屈服面交点在静水压力轴上的投影作为映射中点，并修改剪胀函数。GAO等[86]将描述超固结特性的参数引入修正剑桥模型的剪胀方程中，以模拟在“干侧”发生的剪胀和在“湿侧”发生的剪缩响应。同时，由于临界状态边界面与当前加载面重合，该模型能自然满足经典临界状态理论的要求。

自然沉积的黏土往往还呈现出组构各向异性的特征。各向异性一般可分为两类：固有各向异性和诱导各向异性[87-89]，前者往往是由初始的沉积过程形成，后者则源自于应力状态变化。试验结果[90-91]表明：经过K0固结的土体，其屈服应力轨迹点连线在应力空间中是一个倾斜的椭圆。因此，许多本构模型[84-85,92-94]通过引入旋转硬化机制来描述屈服面倾斜特性，以反映土体应力各向异性程度，并采用相应的旋转硬化法则描述在剪切过程中各向异性的演化。因此，如何确定临界状态旋转角的取值，合理选择旋转硬化法则，并满足临界状态线唯一性和防止屈服面过度旋转是当前研究的热点与重点。DAFALIAS和TAIEBAT[95-96]认为由于临界状态土体具有组构各向异性，因此临界状态旋转角不为零，同时也对临界状态旋转角为零的情况进行了讨论。由于无法通过物理实验证实，该观点尚存争议。COOMBS[97]在自由能中采用了非线性弹性，提出了各向异性超塑性（Hyper-plasticity）模型，该模型在临界状态时的屈服面并非倾斜，旋转角为零。CHEN和YANG[84]从热力学定理出发，认为达到临界状态时由于自由能不能增加，偏背应力必须为零，屈服面转回与静水压力轴平行的位置，而这并不代表土体在临界状态就表现为各向同性，仅是一种模拟上的妥协。微观试验[98-99]和离散元结果[8,13,32,40]表明土体在临界状态呈现出高度的各向异性。然而，基于热力学的本构模型要求旋转角在临界状态时其值为零。CHEN和YANG[84]认为其根源在于传统的旋转硬化机制中旋转角既反映了应力各向异性引起的屈服面旋转，同时也描述了组构各向异性。他们认为若通过引入新的组构张量将旋转角的双重功能解耦，则可以同时满足临界状态旋转角为零的热力学要求和土体在临界状态的高度各向异性特征，相关工作正在开展中。FUENTES等[100]指出如果采用屈服面的旋转角描述土体的固有各向异性，则不能描述与加载方向和沉积平面相关的最大剪切刚度。

除上述对各向异性黏土在单调荷载作用下的研究，许多学者开展了各向异性黏土在循环荷载作用下的研究。WICHTMANN和TRIANTAFYLLIDIS[101]开展了对高岭土进行的循环试验，观察到在不排水循环作用下，达到破坏或者较大应变所需的圈数与沉积平面的方向有关，且q-p循环应力路径是倾斜的。FUENTES等[100]认为这是由于黏土的固有各向异性所致，且所采用的旋转硬化模型无法合理地反映沉积方向对孔压累积、应力路径等的影响。其原因主要是传统的旋转硬化模型[102]不能合理地反映组构方向与加载方向的共同作用，导致有效应力（孔压的累积）模拟不准确。SHI等[103-104]在边界面框架下提出了混合流动法则，以实现黏土在不排水循环中产生的“蝶形”应力路径和合理的孔压累计。

在现有的黏土本构模型中，组构各向异性通常通过屈服面的倾斜角大小考虑。在三轴空间中，该方法被等效为应力比（标量），从而能被纳入多种形式的屈服面（或塑性势面）表达式中。然而，该方法在描述组构方向与加载方向之间的联系时存在一定的局限性。如何使旋转角既能反映应力各向异性（如K0固结）的影响，又能体现由于沉积过程形成的组构各向异性是未来的重点研究方向。

4.2 亚塑性本构模型

近年来，亚塑性理论应用于黏土本构模拟中。在VON WOLFFERSDORFF亚塑性模型[62]的基础上，MAŠÍN[105]提出了首个较为适用的黏土亚塑性模型。MAŠÍ N[106]在亚塑性模型的线性项中引入反映组构各向异性的变量，以模拟黏土小应变刚度的各向异性。与 MAŠÍ N[105]类似，TAFILI和TRIANTAFYLLIDIS[107]和 FUENTES 等[100]在 ISA（Intergranular Strain Anisotropy）亚塑性模型的框架下通过弹性张量来反映沉积平面朝向对土体性质的影响，模拟了循环荷载下黏土固有各向异性对其力学特性的影响。

5 结论与展望

天然土体具有各向异性，对土的强度、剪胀和破坏模式等力学特性具有显著影响。围绕各向异性饱和砂土和黏土，本文介绍了基于弹塑性和亚塑性理论框架土体本构模型的研究进展。

经典临界状态理论为建立能模拟不同密实度的土体本构模型奠定了重要基础。各向异性临界状态理论回答了争论多年的临界状态唯一性问题，为各向异性土体本构模拟提供了理论依据，提升了模型的模拟能力。

针对各向异性砂土，在弹塑性模型中通常将组构各向异性张量与应力张量的第一联合不变量引入屈服函数中，将流动法则、剪胀方程、非共轴响应等自然地与组构各向异性和加载方向联系起来。对于亚塑性模型，组构各向异性变量被引入到特征孔隙率、临界应力面或塑性流动中。

针对各向异性黏土，基于弹塑性理论的本构模型通过屈服面的旋转角大小来表征土体组构各向异性的大小。亚塑性本构模型通过在线性项中引入反映组构各向异性的变量。

尽管众多学者在各向异性土体的组构演化及本构模拟方面取得了显著的进展，但建立模拟更准确、物理意义更明确、普适性更强的本构模型仍面临诸多挑战。对于拟开展的研究，先进实验和数值模拟技术的发展使获得土体在不同应力路径下的组构演化成为可能，这将极大地改善当前本构模拟中采用过多假设的缺陷，特别是从微观角度，建立功能和耗散的统计规律及微观演化机理，获得土体的屈服和硬化规律。而离散元法则是一个非常好的研究工具，为实现上述目标成为可能。