岩土本构综述.docx

第 卷第 期 (小5号宋体) 岩 土 力 学 Vol. No.2008年 月 Rock and Soil Mechanics . 2008文章编号（黑体加粗）：10007598(2003) 02030403（编号用Times New Roman）饱和土本构模型研究进展摘 要：自20世纪50年代以来，随着计算机技术的发展，许多能够描述饱和土体复杂力学行为的本构模型相继被提出来，但由于模型数量较多，很多模型较为复杂，因此不被工程师们所接受。综述近60年来饱和土体静力本构和动力本构的发展情况，对每种模型进行简单的介绍，以求尽可能多的囊括近年来较为成熟的各类模型，便于工程师与科研工作者对这些模型有所了解，并能在工程中进一步完善和应用。关 键 词：中图分类号：TU 443（Times New Roman） 文献标识码：AAdvance in research on constitutive model of saturated soilAbstract: since 1950s, with the development of computer science, many constitutive models were proposed to describe the complicated nature of saturated soil. However, the number of the new model is too large and many of them are not accepted by engineers. We review the development of saturated soil constitutive and soil dynamics constitutive in nearly 60 years, and introduce as many relevant maturity models briefly as possibly in order to make engineers and scientists know about these models and utilize them in real projects.Key words:1. 引言土作为一个自然形成的天然材料，具有复杂的物理力学性质，普遍认为用统一的土的本构模型完全模拟土的物理力学性质是十分困难的1，现有的模型普遍都具有局限性。土体依据颗粒大小，矿物组成等物理性质分为有粘性土与无粘性土，而两类土在力学性质上有很大的不同，尤其是其作为多孔介质材料时，与水发生相互作用，其表现出来的力学性质更是相差甚远。对于同种土不同的含水量也会影响土的力学性质。因此多年来，为了能够较细致的描述土的力学性质，人们一直在针对不同的土给出不同的力学模型，而研究对象也逐渐从饱和土到非饱和土过度。为了适应与更广泛的工程应用，统一的力学模型也是必不可少的。人们运用连续体力学，多孔介质材料力学与混合物理论，给出了土体运动和变形所要满足的各类平衡条件，为了进一步对土体的具体的力学特性进行描述，还需要建立土体的本构方程。对于材料的本构关系的论述最早可追溯到胡克定律，而摩擦型材料需要在线性广义胡克定律的条件下，给出描述摩擦型材料力学特性的莫尔库仑准则。人们最初将土视作为摩擦型材料，因此莫尔库仑模型在很长一段时间被应用到各类岩土工程问题中，直到现在，人们仍然视莫尔库仑准则为土体的破坏准则。在计算机尚不发达的年代，莫尔库仑型理想弹塑性本构模型作为能够模拟摩擦型材料剪切特性的模型起着主导的作用。随着试验技术的发展和越来越多的高精度试验设备的开发，土体越来越多的特性被人们所了解，比如剪胀性，各向异性，结构相关性以及非饱和土的特性在近几十年受到广泛的关注。计算机的发展使得人们可以使用更为复杂的非线性本构关系来描述原本使用莫尔库仑理想弹塑性模型无法描述的土体力学特性2。但是许多很好的模型并没有在工程中得以应用，在进行有限元分析时存在诸多问题。本文将对过去几十年来较为成熟的饱和土体静动力本构模型的研究状况进行简单介绍，以便于更多的工程师对这些模型有所了解，并将这些模型应用于实际工程中去。2. 土体静力本构模型研究进展土体静力本构模型建立了土体在受到静态荷载作用下应力与应变的关系，对于不同的土体，因其密度，受力状态，排水条件等的不同其表现出的应力应变关系有很大的不同3。因此，往往人们在建立土体的本构关系时，会分析其适用条件。土体静力本构模型因其描述的力学特性与模型采用的假设不同，大致可以分为弹性本构模型，弹塑性本构模型，粘弹塑性本构模型，脆性模型和损伤模型等，不同的模型对土体的性质描述的侧重点不同，在采用具体的模型类型时，需要根据实际工程问题来选取。本文主要介绍相对简单并且发展较为成熟的前两类模型。2.1 非线性弹性本构模型非线性弹性模型的共同特点是材料的弹性模型参数与应力应变水平之间有关，对线弹性的广义胡克定律进行了推广，建立增量型、全量型或者积分型的本构关系。依据基本假设的不同，又可以分为：Cauchy弹性模型、超弹性模型和亚弹性模型，相对于Cauchy弹性和超弹性模型，亚弹性模型表达式更具一般性，能够模拟材料更复杂的行为。不同类型的模型定义的应力应变关系函数也不相同。对于Cauchy弹性本构，需要定义了一个二阶张量函数，亚弹性应力应变关系为四阶张量，而超弹性需要定义一个标量势函数，因此从满足能量守恒原则的角度来说，超弹性本构符合热力学原理4。但是不论是Cauchy弹性，超弹性还是亚弹性，待定的材料参数都很多，而且很多参数目前还不能直接由试验给出或者需要进行大量的试验才能得出，大多数的此类模型都只停留在理论层面或者只能通过数值试验来进行参数反演，而这一过程十分复杂，并不被工程界所推崇，因此人们在处理弹性问题时会假设材料为各项同性体或者横观各向异性体来减少本构参数，产生了许多较为实用且参数容易获得的本构模型。典型的非线性弹性本构主要有双曲线型的E-B和E-模型，K-G模型以及小应变刚度模型等。K-G模型是双线性弹性模型的拓展，模型中考虑了应力水平对模型参数的影响，假设模型参数K,G为应力水平的函数，建立增量应力应变关系。典型的K-G模型有：Domaschuk-Valliappan模型，Naylor.模型和沈珠江模型，模型的不同之处主要是对K,G的应力水平函数的假设上。相对于K-G模型，双曲线型模型建立的是全量型应力应变关系，为了能够适应有限元分析，将其化为增量的形式。该模型最初由Kondner(1963)5提出，后来有Duncan进一步发展，得到著名的邓肯张模型(1970)6，并在工程界得到广泛的应用。由于测量手段的局限性，对于微小变形问题其处理方式在很长一段时间都是采用的线弹性假设，随着电子测量技术的提高，人们可以对微小变形进行精准的测量，并发展出一套非线性弹性本构来模拟，Jardine7最早提出了模量随应变水平变化的模型，其后Puzrin8等人对理论进行了完善，提出了微小应变区域（SSR）的概念。以上这些模型大多都是假设材料为各项同性材料，对于像土这种具有横观各向异性的材料，需要对本构关系进行调整，独立的材料参数将会增加，另外，由于上述模型都假设了主应力主应变方向一致，所以对于土体的一些特性如剪胀性不能够很好的模拟，而对于剪胀性等土体特性的模拟，将在后面弹塑性模型中加以介绍。2.2 弹塑性静力本构模型 如上面所述非线性弹性本构模型不能很好地模拟应力应变方向不一致或者状态演化等这些土体相对于金属材料更加复杂的特性，而弹塑性模型可以很好的模拟这些特性，因此近年来受到研究人员的广泛关注。传统弹塑性理论将应变分为可恢复的弹性应变与不可恢复的塑性应变，通过加卸载试验可以确定这些物理量之间的关系，随着增量塑性理论的发展，复杂的土体特性如结构性、剪胀性、各向异性、状态相关性等等都能够通过对增量塑性模型框架进行修改以得到很好的描述，许多成熟的模型涌现出来。其后内变量理论的发展促使人们逐渐弃用了弹塑性应变假设，而提出了诸多直接控制状态变化的内变量来描述更多的力学特性。 塑性本构框架需要满足一个假设，即坐标轴方向一致假设，以及三个准则，即屈服准则，流动法则，硬化准则。它们分别给出了塑性应变增量产生的条件、方向和大小。诸多的弹塑性模型之间的区别主要是采用了不同具体形式的准则，而这些准则是用来描述不同的具体问题而给出的9。例如，使用总应力分析法时分析土的快剪特性往往使用的屈服准则为Tresca准则，而考虑固结影响时因为要考虑固结围压对剪切强度的影响，使用Mohr-Coulomb准则更为合适。但是许多屈服准则在应力空间中的表达式是分段表达式。这给数值分析带来了很多不便，为了处理使用分段函数带来的数值困难，Drucker与Prager于1952年提出了D-P准则10。其准则在应力空间中是一个与Mohr-Coulomb准则相关的圆锥，虽然这样做可以解决屈服面不光滑所带来的数值麻烦，但却不能给出不同方向强度不同的性质，所以为了更好的模拟土体的强度特性， Matsuoka和Nakai(1974)11以及Lade等(1977)12相继提出了基于Mohr-Coulomb准则又能很好的满足光滑性原则的SMP破坏准则。关于流动法则，金属由于其塑性应变的增量方向与应力方向一致，因此往往采用关联流动法则，由于土体塑性应变发展方向并不一定与应力方向一致，采用关联流动法则并不合适。土体剪切屈服面是一个倾斜的曲面，使用关联流动法则会产生过大的剪胀，因此需要采用非关联流动法则来给出土体塑性应变的发展方向。使用非关联流动法则时，需要假设土体存在一个塑性势，塑性应变的发展方向由塑性势面的外法线方向给出，通常认为塑性势与剪胀角有关，这样就可以有效的控制剪胀量的大小而不必受制于屈服面形状的约束。然而，土体是否真的存在塑性势仍然是一个值得争议的问题，对于各向同性不可压缩材料，Hunter(1976)证明了确实存在这样一个势函数，但是对于土体这样的可压缩材料并没有给出一个合理的证明，而对于高维问题，Vardoulakis 和Sulem(1995)也给出了相关的讨论。现在普遍认为高维情况下，仅用一个塑性势函数并不合适。传统的适用于金属材料的应力应变分析方法不能很好地模拟土体排水与不排水特性，随着有效应力分析法的广泛应用，弹塑性本构模型使用有效应力表示。由于决定土体力学性质的状态量不仅仅有应力应变，还有土体的物理指标，比如密度、颗粒形状、沉积情况等等。20世纪50年代，针对状态变化影响土体的力学特性，人们展开了临界状态模型的研究，将土体的压缩特性与剪切特性在模型中统一起来，建立了一系列能够很好描述岩土材料力学特性的本构模型。 从硬化参量选取的角度可以将硬化规律分为应变硬化与功硬化两类，三轴状态下应变硬化更容易让人理解，因此，在描述土体的三轴状态硬化软化特性时，往往采用塑性体积应变硬化和剪切应变硬化。Drucker(1957)13等提出帽子屈服面的概念，将屈服与体积应变联系起来，Dimaggio和Sandler(1976)基于Drucker等人的研究提出帽盖模型，使模型能够描述体积硬化，体积硬化规律往往由固结曲线给出，剪切硬化规律则由剪应变等值面硬化规律给出。从加载路径的角度又可以将硬化规律分为等向硬化、运动硬化、混合硬化以及旋转硬化四类。对于土体排水条件下的各向同性等体压缩，采用等向硬化模型即可，而当考虑循环加卸载Bauschinger效应时，需要引入运动硬化来保证应力应变滞回的封闭性。对于K0条件下的土体压缩问题，土体会表现出拉压各向异性，屈服面在应力空间中会发生旋转，因此需要采用旋转硬化规律。剑桥大学的Roscoe最早就正常固结土和超固结土剪胀性问题提出了著名的剑桥模型。在该模型框架中，引入比体积这一土体的状态参量，假定存在状态边界面。屈服面在临界状态的左面和右面各自给出土体不同的力学性质，当进行排水分析时，屈服面的演化由先期固结压力来控制。由于剑桥模型在偏应力等于0处存在奇异点，对数值和理论分析带来不便，因此Roscoe等人有发展了剑桥模型，提出了基于椭圆形屈服面的修正剑桥模型，其后的几十年里，人们不断地改造和发展了临界状态弹塑性模型，使之能够描述更多的土的力学行为。Tanaka(1986)14采用Hvorslev(1937)15面为干面，并采用非关联流动法则，以解决剑桥模型与修正剑桥模型计算所引起的过高的剪胀率以及过高的峰值强度。土的各向异性可以分为土体固有的各向异性（初始各向异性）和当土体受力时, 应力不等引起的各向异性（诱发的各向异性），后者对土体力学性质影响更大16。Ohta(1976)17，Kavvadas(1982)18以及Sekiguchi(1977)19等人将模型的屈服中心移至K0线上，使得屈服面在应力空间中发生了偏转，以使模型能够模拟偏压固结时土体表现出来的诱发各向异性的特性。而后Mouratidis(1983)20提出了MELANIE模型，能够分析非关联流动K0固结问题，由于沿K0状态方向剪切时，使用剑桥模型计算的体积应变偏小，而使用关口太田模型计算的体积应变偏大，所以孙德安和姚仰平(2000)21提出了介于上述两者之间并能考虑初始各向异性的K0硬化模型。为了能模拟屈服面内应力应变非线性。Hashicuchi(1977)22等人提出次加载面模型，Mroz(1978)23等提出出多重屈服面模型，Dafalias(1982)24等提出边界面模型，Pande(1982)25等人提出广义塑性模型，这些模型共同特点是认为边界面内存在控制应力应变非线性发展的屈服面，该屈服面与边界面存在几何相似性，Mroz的多屈服面模型需要记录大量的屈服历史，而其他模型只需要记录内面与外面的相对位置即可。Perzyna等(1963)25提出能够较好的描述应力应变时间关系的超应力模型来描述软土的次固结性质。Hsieh等(1990)26在剑桥模型中加入了时间相关项，使得其可以允许应变延迟。Rouainia和Wood(2000)27将时间影响项引入到边界面模型的理论框架。土体的正常固结与次固结性质往往与其结构性的改变有关，所以许多研究者就土体结构性展开了一些列的研究。Liu和Carter(2002)28对修正剑桥模型进行了改造，是剑桥模型能够描述土的结构对土的力学性质的影响。通过对原状土力学性质的研究，人们认为土体的变形过程看做是由原状土经损伤向扰动土逐渐转化的过程1，因此可以采用损伤力学建立土体的弹塑性损伤模型来描述这一转化过程。对于土的结构性本构模型归纳来说有三种研究途径：一是微观形态学的研究，二是固体力学研究方法，三是土力学的研究方法，对于微观形态学主要考察土体变形过程中的微观结构的形态变化，利用非线性动力学方法来研究土体微结构参数与宏观连续介质力学参数之间的联系，建立微观力位移与宏观应力应变的关系。利用损伤力学方法来研究土的结构性，由于是基于宏观力学方法，所以相对更容易被工程界所接受，沈珠江和章为民(1988)通过定义理想原状土和完全损伤土，而真实土体介于二者之间，给出了一套适合于土体的损伤理论。Desai(1992)29在本构模型中引入了扰动状态的概念。Armaleh(1994)30等进一步发展了Desai的模型，并给出了相应的试验方法。随着组构张量概念的提出与发展，组构张量逐渐替代孔隙率而成为新的描述岩土材料密实状态以及土体结构的物理量，对于结构对土体力学性质的影响近年来许多研究者做了大量的工作，本文不做具体介绍，具体请参阅31-36。增量型本构关系（尤其是非线性弹性本构）虽然能够很好的给出的应力应变关系，但很多情况下不能满足热力学基本定律，违背基本的物理守恒，因此，很多本构假设需要给出其是否满足热力学定律的验证过程。Ziegler最早将热力学方法引入到固体应力应变分析中，提出了Ziegler正交假定，并指出大多数材料满足Ziegler正交假定而不满足Druker公设。Houlsby(1981)最先将热力学方法引入到土体本构理论中，给出了修正剑桥模型对应的自由能函数和耗散函数形式，提出耗散应力空间中屈服面的概念。其后Collins和Houlsby(1997)进一步给出了使用热力学方法构造率无关超塑性本构模型的基本框架。Houlsby 和Puzrin(2002,2003,2006)又给出了率相关超塑性本构框架，对于剪胀性，应力诱发各向异性等这些土体特有的力学特性，Collins和Muhunthan (2003)应用热力学方法分析了在耗散应力空间中描述剪胀性的剪胀角和描述各向异性的各向异性角的关系。3. 土体动力本构模型研究进展饱和土体在动荷载作用下的应力应变关系是极其复杂的，它取决于土性条件、荷载条件、排水条件等等，与土体的静力模型一样，想要建立一个能够满足所有条件的动力本构模型是很难做到的，应当根据实际工程问题，选取主要的影响因素，建立能够反映实际问题的本构模型。土体在动荷载作用下的变形包括弹性变形与塑性变形，对于循环应力作用是其应力应变关系具有非线性、滞回性和棘轮效应等特点。在早期进行的小变形研究时，人们会假设土为近似弹性体或粘弹性体，研究剪切模量与阻尼比的变化规律，而对于较大的变形，土体结构会发生破坏，孔压上升致使强度丧失发生液化。对饱和土体进行液化分析时，往往需要将土体视为弹塑性体。本文将重点对动力粘弹性模型与弹塑性动力本构模型进行介绍。3.1 动力粘弹性本构模型 土体的动应力-动应变曲线主要由骨干曲线和滞回曲线构成，骨干曲线给出的是动应力和动应变的非线性，而滞回曲线给出了动应力和动应变的滞后性，正如前面所说，非线性和滞后性是主要的两个土体动力特性，因此人们在很长一段时间都在研究土体的动力粘弹性本构，并给出了很多实用的模型。 土体动力粘弹性本构模型的建立最早应追溯到1968年Seed在分析土体动力非线性性质所使用的等价线性模型，模型并没有给出骨干曲线和滞回曲线的具体表达式，而是将动力的非线性归结为模量随剪应变的衰减，将滞后性归结为阻尼比随剪应变的升高。随着粘弹性模型的发展，卸载再加载的曼辛效应逐渐引起研究者的注意，滞回效应不再是简单的使用阻尼比来控制，而是更加关心滞回曲线的形态。比如Hardin-Drnevich模型、Ramberg-Osgood模型以及其他的组合曲线模型，一般称这类模型为Masing类模型。Pyke(1979)对Masing准则做了修正，使其能够适用于非等幅作用情况下的问题。Ishihara等(1985)对Masing准则进行了修正，使其能够描述更大应变范围的动力性质。一般的等价粘弹性模型不能计算永久变形，因此Martin等人根据不排水孔压和排水体应变之间的关系，给出能计算永久变形的模型，除上述的基于试验数据拟合经验公式的模型外，伊万等人提出元件模型，即使用弹簧阻尼塑性元件的串并联来模拟应力应变的非线性，而土工试验作为一类特殊的边值问题来考虑。郑大同和王惠昌(1983)，引入了软化系数，使Iwan模型既能反映硬化又能反映软化，之后李小军和廖振鹏(1989) 以1wan并联模型为基础提出了一种具有粘性、弹性及塑性三种效应的应力应变关系物理模型粘一弹一塑性三元素模型。动力粘弹性本构模型由于其相对简单，有大量试验数据支撑并其模型参数可由试验获得，在工程界得到广泛的认可和应用，但其仍然有很多缺陷，比如不能考虑应力路径影响，不能考虑密度变化影响以及应力诱发各相异性，其永久变形分析需要配合孔压分析，人们为此对粘弹性本构做了很多的修正，使得模型非常复杂，不利于推广。因此，动力弹塑性模型在这些方面具有优势。3.2 弹塑性动力本构模型 动力弹塑性模型相对静力弹塑性模型需要考虑动力作用下的滞回效应以及棘轮效应。传统塑性力学中往往使用运动硬化来描述应力应变滞回，即屈服中心会在应力空间中随着加卸载过程发生变化，进一步考虑运动硬化与等向硬化相结合的混合硬化37-39。随着临界状态土力学的发展，人们将临界状态模型进行改造使其能够模拟岩土材料的动力特性，这样弹塑性模型便不区分静力模型还是动力模型，动力特性可以通过作用于特殊动应力路径给出的结果。Mroz(1967)首次提出塑性模量场，其基本概念为在应力空间中存在一个边界面，边界面以内存在一些列相互嵌套但不相交的嵌套面，每个嵌套面对应一个硬化模量，这样在应力空间中便存在了一个模量场，可以描述加卸载的非线性应力应变关系以及滞回效应，之后，Mroz等人进一步发展了该模型，使其能够适用于描述排水条件下的砂土振动压密的现象，同时可以考虑扭剪和体变。Prevost等40, 41研究了粘性土在排水和不排水条件下的应力应变强度特性，并建立了能够模拟该特性的多屈服面模型。Prevost等，Zienkiewicz等42进一步发展了多屈服面模型。为了使应力应变关系光滑性更好， Mroz提出了无穷屈服面模型。多屈服面模型和无穷屈服面模型虽然在模拟土体滞回应力应变非线性具有更好的灵活性和普遍性，但是需要对所有的屈服面的位置、尺寸以及塑性模量进行记忆，这对于进行数值计算是非常不利的，会消耗大量的内存。因此Dafalias和Krieg分别提出了两面模型，即边界面内部的屈服面是与边界面几何相似，该屈服面遵循一定的规则移动，屈服面上应力点处的硬化模量可由其对应在边界面上的共轭点处的硬化模量和两点的相对位置插值得到，而边界面上的塑性模量与外法线方向由塑性硬化规律和流动法则给出。Poorooshasb等(1985)基于边界面模型框架，给出了能够模拟三维情况单调加载和循环加载的本构模型，其后有进一步发展为能够模拟复杂应力路径的本构模型。Desai(1985)和Hirai(1987)修正了两面模型使其能够描述应力诱发各向异性。Tabaa和wood(1989)对修正的剑桥模型进行了扩展，建立了一个用于描述粘土在循环荷载作用下滞回特性的双面动力硬化模型。Cmuch和wolf(1994)建立了适用于粘土、粉土和砂土的统一的三维临界状态边界面塑性模型，Borja和Alnies(1994)修正三维边界面模型使其能够描述粘土在多轴应力反向条件下的应力应变特性， Dafalias和Manzari(2006)43研究了粘土的旋转硬化特性，并给出了能够更加准确的描述粘土各向异性的边界面模型SANCLAY。Wang等(2013)44修正了边界面模型，使其能更好的描述砂土不排水相变时的剪胀性，模型计算结果与实验结果能够很好得吻合。由Hashiguchi(1985)提出的次加载面模型，可以很好地克服两面模型在数值计算是出现屈服面超出边界面的情况。通过相似中心的迁移和相似比的改变，给出次加载面与正常屈服面的关系及二者的运动规。HashiguchiWhittle(1987)基于边界面模型提出了MIT模型，该模型能够描述土体的各向异性，循环特性以及不排水情况下的应变软化特性等等。而其后发展的MIT-E3模型采用了边界面理论并可以准确的反映小应变情况下的闭合应力应变滞回,能够很好地模拟卸载再加载的塑性应变积累.4. 展望虽然土体本构模型发展已有数十年的历史，由于土体力学性质过于复杂，仍然有很多问题不能很好地解决，人们一方面不断的完善过去几十年来建立的本构模型，是它们能够描述更多的力学性质，另一方面也在尝试着从不同的角度来分析土体的力学行为。近十年来兴起的微细观力学以及土体结构化的研究越来越受到科研人员的重视，人们从微观角度对饱和土体力学行为进行定性的解释，并结合以DEM（离散单元法）进行定量分析，建立微细观与宏观力学行为的联系。虽然微细观机理的研究成为趋势，但是工程上仍然习惯于使用宏观本构模型进行工程分析，因此在完善宏观本构理论和机理研究的同时，应当大力地将这些成果在工程中进行推广。References: 1.张建民与谢定义, 饱和砂土动本构理论研究进展. 力学进展, 1994. 2. 2.郑颖人, 沈珠江与龚晓南, 岩土塑性力学原理. 2002: 中國建筑工业出版社. 3.章根德, 土的本构模型及其工程应用. 2003: 科学出版社. 4.李广信, 高等土力学. 2004: 清华大学出版社有限公司. 5.Kondner, R.L., Hyperbolic stress-strain response: cohesive soils. Journal of the soil mechanics and foundations division, ASCE, 1963. 89(1): p. 115-143. 6.Duncan, J.M. and C. Chang, Nonlinear analysis of stress and strain in soils. Journal of the Soil Mechanics and Foundations Division, 1970. 96(5): p. 1629-1653. 7.Jardine, J.R., et al., Studies of the influence of nonlinear stress-strain characteristics in soil-structurs interaction. Geotechnique, 1986. 3(36): p. 377-396. 8.Zdravkovi, O.M.P.A., Finite Element Analysis in Geotechnical Engineering - Theory. 9.郑颖人, 岩土弹塑性力学. 2010, 中国 北京: 中国建筑工业出版社.10.Drucker, D.C. and W. Prager, Soil mechanics and plastic analysis of limit design. Quarterly of Applied Mathematics, 1952(10): p. 157-164.11.Matsuoka, H. and T. Nakai. Stress-deformation and strength characteristics of soil under three different principal stresses. 1974.12.Lade, P.V., Elasto-plastic stress-strain theory for cohesionless soil with curved yield surfaces. International Journal of Solids and Structures, 1977. 13(11): p. 1019-1035.13.Drucker, D.C., R.E. Gibson and D.J. Henkel, Soil mechanics and work hardening theories of plasticity. Trans.ASCE, 1957. 122: p. 338-346.14.Tanake, T., M. Yasunaka and Y.S. Tani, Seismic response and liquefaction of embankments-numerical solution and shaking table tests, in Int.Symp.Num.Mod.in Geomech. 1986: Ghent. p. 679-686.15.Hvorslev, M.J., Uber die Festigkeitseigenschaften gestorter bindiger Boden. 1937.16.殷宗泽与徐志伟, 土体的各向异性及近似模拟. 岩土工程学报, 2002. 24(5): 第547-551页.17.Ohta, H. and C.P. Wroth. Anisotropy and stress reorientation in clay under load. 1976.18.Kavvadas, M., Non-linear consolidation around driven piles in clays. 1982, Massachusetts Institute of Technology.19.Sekiguchi, H. and H. Ohta, INDUCED ANISDTROPY AND TIME DEPENDEHCY IN CLAYS. 1977.20.Mouratidis, A. and J. Magnan, An anisotropic elastoplastic model with increase in strength through plastic deformation for the design calculation of structures on compressible soils. 1983.21.孙德安与姚仰平, 初始应力各向异性土的弹塑性模型. 岩土力学, 2000. 21(3): 第222-226页.22.Hashiguchi, K. An expression of anisotropy in a plastic constitutive equation of soils. 1977.23.Mrz, Z., V.A. Norris and O. Zienkiewicz, An anisotropic hardening model for soils and its application to cyclic loading. International Journal for Numerical and Analytical Methods in Geomechanics, 1978. 2(3): p. 203-221.24.Dafalias, Y.F., On cyclic and anisotropic plasticity: i) a general model including material behavior under stress reversals: ii) anistropic hardening for initially orthotropic materials. 1975, University of California, Berkeley.25.Pande, G.N. and S. Pietruszczak. Reflecting surface model for soils. 1982.26.Hsieh, H.S., E. Kavazanjian Jr and R.I. Borja, Double-yield-surface Cam-clay plasticity model. I: theory. Journal of Geotechnical Engineering, 1990. 116(9): p. 1381-1401.27.Rouainia, M., A kinematic hardening constitutive model for natural clays with loss of structure. Geotechnique, 2000. 50(2): p. 153-164.28.Liu, M.D. and J.P. Carter, A structured Cam Clay model. Canadian Geotechnical Journal, 2002. 39(6): p. 1313-1332.29.Desai, C.S. and Y. Ma, Modelling of joints and interfaces using the disturbedstate concept. International Journal for Numerical and Analytical Methods in Geomechanics, 1992. 16(9): p. 623-653.30.Armaleh, S.H. and C.S. Desai, Modelling and Testing of a Cohcsionless Material Using the Disturbed State Concept. Journal of the Mechanical Behavior of Materials, 1994. 5(3): p. 279-296.31.Pietruszczak, S. and S. Krucinski, Description of anisotropic response of clays using a tensorial measure of structural disorder. Mechanics of Materials, 1989. 8(2): p. 237-249.32.Kachanov, M. and I. Sevostianov, On quantitative characterization of microstructures and effective properties. International Journal of Solids and Structures, 2005. 42(2): p. 309-336.33.Mehrabadi, M.M., S. Nemat Nasser and M. Oda, On statistical description of stress and fabric in granular materials. International Journal for Numerical and Analytical Methods in Geomechanics, 1982. 6(1): p. 95-108.34.Bathurst, R.J. and L. Rothenburg, Observations on stress-force-fabric relationships in idealized granular materials. Mechanics of Materials, 1990. 9(1): p. 65-80.35.Ken-Ichi, K., D