基于微观力学分析的散粒体静力液化本构模型

刘 洋，汪成林，闫鸿翔

（1.北京科技大学 土木与环境工程学院土木系，北京 100083；2.中国公路工程咨询集团有限公司，北京 100097）

基于微观力学分析的散粒体静力液化本构模型

刘 洋1，汪成林1，闫鸿翔2

（1.北京科技大学 土木与环境工程学院土木系，北京 100083；2.中国公路工程咨询集团有限公司，北京 100097）

针对散粒体的静力液化特性，在散粒体颗粒运动微观力学分析的基础上，在临界状态土力学框架内建立了一个简单的静力液化弹塑性本构模型。模型屈服函数和和硬化规律根据Chang提出的砂土微观力学模型，通过积分粒间接触力为宏观的应力不变量而建立，考虑了与材料状态相关的剪胀性和初始密度对散粒体应力 应变关系的影响，并采用非相关联的流动法则。模型参数简单且有较明确物理意义，可以通过室内三轴试验确定。模型的数值结果与Toyoura砂以及砂 粉混合土的三轴不排水剪切试验的应力 应变曲线和应力路径吻合较好。

散粒体；静力液化；本构模型；数值模拟

砂土等散粒体的液化包括静荷载引起的静力液化和地震、爆炸以及机械振动等动力荷载引起的循环液化。对于循环液化，学者们对其进行了大量的研究［1］；对于静力液化，近几十年间却没有得到足够的关注和认识。对地基、边坡等土工构筑物来说，静力液化是一种常见的破坏形式，对土工问题进行静力液化稳定性分析时，不管采取何种数值方法，首要问题是建立合理的静力液化应力－应变关系。

散粒体（如砂土）的应力 应变模型主要有两类：基于连续介质力学的弹塑性本构模型和基于散粒介力学分析的微观力学模型。经典弹塑性力学模型将砂土看成连续介质，如Desai等［2］提出的帽子模型、Prevost等［3］提出的多重屈服面塑性模型以及Dafalias等［5］提出的边界面模型等等。另外一些学者将砂土看成颗粒集合体，基于接触力学和均一化方法提出砂土的微观力学模型，早期的研究主要集中在砂土弹性力学行为，如Chang等［5］的工作。对弹塑性应力－应变关系的研究要复杂一些，如Jenkins等［6］、Chang等［7］、Yin 等［8］、Nicot等［9－10］、Misra等［11］、Lai等［12］、Zhang等［13］、Tran 等［14］、Shen等［15］提出的微观力学模型。对于散粒体的静力液化问题，也有学者提出了一些用以模拟各向同性不排水条件下砂土的静力液化的模型（如Boukpeti等［16－17］）。

本文针对散粒体的静力液化特性，借鉴微观力学本构关系一些成果，在临界状态土力学框架内建立散粒体的静力液化本构模型。文章首先简要阐述了静力液化的基本原理，接着提出了适用于散粒体的静力液化本构模型，并对模型参数及其一般力学响应进行了分析，最后根据已有的三轴试验结果对提出的模型进行了评价。

1 静力液化

关于静力液化，太沙基早在1948年首先用“自然液化”描述了非常松散的砂土在微小扰动下突然失去强度而象粘滞性流体一样的流动现象，这里的“自然液化”即为后来 Castro等［18］和 Casagrande［19］在讨论滑坡灾害时提出的“静力液化”概念。它们描述了在静荷载条件下砂土的强度降至很低，不能再继续承受剪切作用，可以如液体一样流动的特性，同时也指出只有非常松散的砂土才能发生静力液化。

从概念上静力液化可以定义为：在单调加载过程中偏应力－应变曲线出现明显的应变软化现象，随着偏应力在峰值后的急剧降低而接近零值，砂土表现出类似于流体的特征，但由于砂土未承受动荷载作用，为了与通常的振动液化区别，故称之为“静力液化”。

静力液化的原理如图1所示，τst是砂土的初始应力状态，φcv是临界状态摩擦角，曲线①、②分别是单调荷载和循环荷载下的应力路径。无论是单调荷载还是由振动产生的循环荷载，由它们产生的瞬态荷载都会使砂土产生液化。图中可见，应力一旦超过破坏面，剪切强度就会快速下降，至到稳定状态即残余强度sus处。虽然单调荷载和循环荷载的应力路径不同，但它们都具有相同的最终状态。

图1 饱和散粒体单调和循环加载响应

2 散粒体颗粒运动微观力学分析

2.1 颗粒间弹性力学行为

定义散粒体颗粒间的接触法向为垂直于接触面的法向矢量，在每个接触面上建立局部坐标如图2所示。颗粒接触法向与切向接触刚度分别为kαn和kαr，粒间接触力可定义为

图2 粒间接触处的局部坐标系

式中，n、s、t是局部坐标系中3个正交的单位法向量。设粒间接触服从修正的Hertz-Mindlin接触本构关系［5］：

2.2 颗粒间塑性力学行为

颗粒接触处存在着向上或者向下的塑性滑动，这种塑性运动产生体积改变，在散粒体中称之为剪胀／剪缩行为，粒间剪胀可用式（4）描述。

3 散粒体静力液化宏观本构方程

基于上述颗粒运动微观力学分析，假设散粒体均匀，并在各向同性应力状态下积分式（4）、（6）、（7），即得到散粒体宏观剪胀方程、屈服函数和硬化规律。

3.1 控制方程

模型的数学表达式以应力和应变不变量的形式给出，应力不变量采用p′和q，它们分别与第一应力不变量I1和第二偏应力不变量J2相关，式中p′和I1′代表有效应力。p′和q以应力张量和偏应力张量表示为

式中，偏应力张量sij＝σii－1／3σkkδij，与其相应的应变增量以应变张量和偏应变张量表示为

其中，偏应变张量eij＝εii－1／3εkkδij。假设材料各向同性，根据弹塑性理论材料变形分为弹性和塑性两部分，即：

这里K0和是两个材料常数，K0是参考体积模量；pref是参考压力，一般取1个大气压；为指数。

上述弹性行为只发生在屈服面内，屈服准则采用摩尔－库伦类准则

该屈服面是在Chang等［14］提出的砂土微观力学模型的基础上，将微观的接触力积分为宏观的应力和而来。同理，根据式（7）硬化规律，采用式（19）所示双曲线形式。

式中，Mp＝6sinφp／（3－sinφp），从物理意义上讲，φp的值应等于库伦摩擦角（φμ），但为了考虑散粒体密实度的影响，消除各向同性假设带来的误差，将φp看成是密实度函数，即

式中：m是一个与散粒体类型（特别是颗粒形状）有关的一个正指数；φcs是临界状态摩擦角。对于密砂，峰值摩擦角φp大于临界状态摩擦角φcs，在剪应力作用下砂土剪胀，颗粒的咬合和摩擦降低，产生应变软化行为。

这里ptm是参考大气压力。式（21）中的3个参数确定了e-p空间临界状态线（CSL）的形状，即：eerf（CSL在e轴的截距），λ（CSL线的斜率），以及ξ（CSL的曲率）。

塑性刚度kp0与体积模量K的关系用式（22）关联。

显然，本文提出的模型采用了非相关联的流动法则。

2.2 应力 应变关系

根据经典弹塑性理论，总应变增量分成两部分，即弹性部分和塑性部分。

4 本构方程的数值积分

在积分计算前，首先假设材料响应为弹性响应，通过弹性本构关系获得一个测试应力状态，即：

式中，T表示应力状态为测试应力，数值大小为当前应力状态加上微小变形后所产生的弹性应力。如果测试应力状态点位于屈服面所定义的弹性区域内，则有效应力分量的更新值等于测试应力。

如果测试应力状态点位于当前屈服面外，即F（p′T，qT）＞0，则对测试应力进行修正计算，此时积分步计算分为两个子步：求解新的应力状态点，使其满足屈服条件；更新屈服面位置，进行更新计算。积分步如图3所示，图中F代表屈服面，G代表塑性势面。

图3 本构方程积分步的两个子步

新的应力状态是在弹性变形和塑性变形共同作用下达到的，如图3（b）所示，因此需要由弹性法则和塑性应变增量共同来确定。

结合流动法则可将新的应力状态表示为式（35）和（36），点 （p′NqN）位于屈服面上。

式（38）的求解采用牛顿迭代法。当计算结果为负值时，塑性因子将取为零，相当于没有发生塑性变形。应力不变量的更新值可由式（35）、（36）计算获得，有效应力张量为

在上述计算完成后，还需采用硬化规律计算屈服面的更新位置，同时对硬化参数进行更新。

式（40）中Mp的计算要采用当前的孔隙比。

5 模型参数分析及力学响应

5.1 模型参数

模型参数总计10个如表1所示，即：3个弹性常数K0、μ（或G）和n；4个塑性参数ζp、D、m 和φcs（φcs也是临界状态参数）；3个临界状态参数eref、λ、ξ。

表1 模型参数表

除了临界状态参数（包括φcs），其余6个参数均可以通过各向同性压缩和三轴剪切试验获得，限于篇幅，这里不再详述参数确定程序。

5.2 模型力学响应

图4是本构模型在不排水三轴应力条件下的力学响应，模型参数取值如表1，图4（a）是不同初始孔隙比试样在同一固结围压（200 k Pa）下的应力 应变曲线，图4（b）是对应的应力路径；图4（c）、（d）是同一初始孔隙比试样在不同固结围压下的应力 应变关系曲线和应力路径。

从图中可以看出，本文提出的本构模型能够反映密实状态（或初始孔隙比）对散粒体力学行为的影响，如图4（a）、（b）所示：在同一固结围压下不同初始孔隙比导致了散粒体不同的剪胀、剪缩行为。不排水条件下，随着剪应变的发展，对于初始孔隙比＝0.7、0.725的密实散粒体发生剪切硬化，而初始孔隙比＝0.75、0.8松散散粒体在达到峰值应力后产生了剪切软化，随着超孔隙水压力的增长，有效应力降低，应力路径向原点发展，发生了静力液化。

模型还能考虑散粒体力学响应的围压依赖性，如图4（c）、（d）所示，对初始孔隙比为0.8的散粒体在不同固结围压（200、400、800、1 200 kPa）下表现出不同的力学响应，低固结围压下散粒体更易于发生静力液化。本构方程的数值结果显示，本文提出的模型可以模拟散粒体如砂土等的不排水剪切力学特性。

6 模型的试验验证

为了进一步验证模型，首先选Toyoura砂三轴不排水试验结果进行比较，然后对砂 粉混合土的三轴不排水试验结果进行模拟，验证本文提出的模型应用于一般散粒体的可行性。

6.1 Toyoura砂

图4 不同初始孔隙比和固结围压的模型响应

Verdugo等［20］对Toyoura砂进行了一系列三轴不排水剪切试验，试验固结围压为100、1 000、2 000、3 000 k Pa。本文对初始孔隙比为0.833的Toyoura砂在不同初始围压下的不排水响应进行模拟，图5是本文模型的预测结果与试验数据的对比，其中连续曲线为数值计算结果，模型参数见表2。

从图5中可见，模型模拟结果和试验结果基本吻合，虽然在应力－应变曲线峰值后的一小段区间内应力水平要略小于试验值，以及应力路径的中部稍有偏差，但从整体趋势看，本文提出的模型较好地描述了Toyoura砂的不排水三轴力学特性。

表2 模型参数表

图5 Toyoura砂试验与模型模拟结果

6.2 砂 粉混合土

Yang［21］对Hokksund砂和中国黄河下三角洲埕北海区粉土混合土进行了一系列不排水三轴试验。试验在3个围压5、100、150 k Pa下进行，除了部分试样在低围压很快液化，其余试验压缩至轴向应变20%左右达到稳态。限于篇幅，我们选取粉粒含量为5%、30%两组试验结果来评价本文提出模型的应用于砂 粉混合散粒体的可行性。模型参数如表2所示，图6是应力 应变关系曲线和应力路径的试验数据和模型预测结果。

图6 砂 粉混合土试验与模型模拟结果

从图中可以看出，对低粉粒含量（5%）的混合土，低围压下更容易产生液化，50、100 k Pa的模拟结果与试验结果较吻合，150 k Pa的试验点较离散，模拟结果稍差；对30%含量的混合土，全部围压下试样均产生了完全的静力液化，有效应力减少至零，模型模拟结果在峰值强度上稍低，但整体上与试验结果吻合较好。

6.3 讨 论

从对Toyoura砂和砂 粉混合土三轴不排水试验的模拟结果看，本文提出的模型可以较好地应用于一般砂土和砂 粉混合散粒体。但需要指出的是，对于砂 粉混合散粒体，基于颗粒形状、砂和粉粒颗粒粒径和相对含量的多少，需要采用不同的模型参数进行预测，从表2中可以看出，对不同粉粒含量的混合土，模型参数虽有差别但比较接近，进一步工作应寻求砂 粉土混合物模型参数与粉粒含量之间的关系，期望能用统一的纯砂和纯粉土参数模拟砂 粉混合散粒体。

此外，本文模型的特点是在散粒体颗粒运动微观力学分析的基础上，模型的屈服函数和剪胀方程是通过积分粒间接触本构方程和运动方程求得，这与一般弹塑性本构模型的建模思路不同。与Boukpeti等［18－19］提出的静力液化本构关系比较，本文的模型参数较少且均具有明确的物理意义，可以通过室内常规三轴压缩和剪切试验结果获取，方便应用。但与Dafalias等［22］提出的本构关系比较，本文的模型考虑了剪切硬化但没有考虑旋转硬化，因此在应用于各向异性砂土力学分析时尚有误差，进一步的工作应考虑建立旋转硬化或者引入组构参数使模型能够应用于各向异性砂土的本构关系的模拟。

7 结 论

针对散粒体的静力液化特性，基于散体微观力学分析在临界状态土力学框架内建立了一个的静力液化本构模型，研究初步得出了以下几点结论：

1）提出的模型采用了临界状态理论，单调荷载下的排水和不排水应力路径在临界状态处终止。模型硬化规律较合理，并采用了非关联的流动法则，考虑了与材料状态相关的剪胀性以及初始孔隙比、固结围压对散粒体应力－应变关系的影响。

2）本构模型的数值结果显示，提出的模型可以较好地描述散粒体如砂土的不排水三轴剪切应力应变关系和应力路径，模型参数简单且物理意义较明确。

3）提出的模型与Toyoura砂和砂 粉混合土三轴不排水试验的结果吻合较好，可以应用于散粒体的静力液化问题分析。

［1］王刚，张建民.地震液化问题研究进展［J］.力学进展，2007，37（4）：575-589.

Wang G，Zhang J M.Recent advances in seismic liquefaction research ［J］.Advances in Mechanics，2007，37（4）：575-589.

［2］Desai C，Siriwardane H.Constitutive laws for engineering materials with emphasis on geologic materials ［M］.Printice-Hall，Eaglewood Cliff，NJ，1984.

［3］Prevost J.A simple plasticity theory for cohesioless soils［J］.Soil Dynamics and Earthquake Engineering，1985，4（1）：9-17.

［4］Dafalias Y，Herrmann L.Bounding surface formulation of soil plasticity ［C］／／Soil Mechanics-transient and Cyclic Loads，Wiley，1982，London：253-282.

［5］Chang C S，Misra A.Initial moduli of particulated mass with frictional contacts［J］.International Journal for Numerical and Analytcal Methods Geomech，1989，13（6）：629-644.

［6］Jenkins J T，Strack O D.Mean-field inelastic behavior of random arrays of identical spheres［J］.Mechanics of Material，1993，16：25-33.

［7］Chang C S，Hicher P Y.An elasto-plastic model for granular materials with microstructural consideration［J］.Intermational Journal of Solids ＆ Structures，2005，42：4258-4277.

［8］Yin Z Y， Chang C S， Hicher P Y，et al.Micromechanical analysis of kinematic hardening in natural clay ［J］.International Journal of Plasticity，2009，25（8）：1413-1435.

［9］Nicot F，Darve F.Basic features of plastic strains：From micro-mechanics to incrementally nonlinear models［J］.International Journal of Plasticity，2007，23：1555-1588.

［10］Nicot F S， Darve F.Failure in rate-independent granular materials as a bifurcation toward a dynamic regime［J］.International Journal of Plasticity，2012，29（1）：136-154.

［11］Misra A，Yang Y.Micromechanical model for cohesive materials based upon pseudo-granular structure ［J］.International Journal of Solids ＆Structures，2010，47：2970-2981.

［12］Lai Y Y，Li S.Strength criterion and elastoplastic constitutive model of frozen silt in generalized plastic mechanics ［J］.International Journal of Plasticity，2010，26（10）：1461-1484.

［13］Zhang W，Zhao C. Micromechanics analysis for unsaturated granular soils［J］.Acta Mechanica Solida Sinica，2011，24（3）：273-281.

［14］Tran T H， Bonnet G. A micromechanics-based approach for derivation of constitutive elastic coefficients of strain-gradient media ［J］.International Journal of Solids ＆Structures，2012，49（5）：783-792.

［15］Shen W Q，Gatmiri B.A micro-macro model for clayey rocks with a plastic compressible porous matrix ［J］.International Journal of Plasticity，2012，36：64-85.

［16］Mro′z Z，Nboukpeti N，Drescher A.Constitutive model for static liquefaction［J］.International Journal of Geomechanics，2003，3（2）：133-144.

［17］Boukpeti N，Drescher A.Triaxial behavior of refined Superior sand model［J］.Computers and Geotechnics，2000，26：65-81.

［18］Castro G，Poulos S，France J，et al.Liquefaction induced by cyclic loading ［R］.National Science Foundation，1982.

［19］Casagrande A.Liquefaction and cyclic mobility of sands，a critical review［D］.Harvard University，1976.

［20］Verdugo R，Ishihara K.The steady state of sandy soils［J］.Soils and Foundations，1996，36（2）：81-91.

［21］Yang L S.Characterization of the properties of sand-silt mixtures ［D］. Norwegian Univ. Sci. Tech.，Trondheim，Norway，2004.

［22］Dafaliasy F，Manzari T M.Simple plasticity sand model accounting for fabric change effects［J］.Journal of Engineering Mechanics，2004，130（6）：622-634.

（编辑 王秀玲）

A Static Liquefaction Constitutive Model for Granular Materials Based on the Micromechanical Analysis

Liu Yang1，Wang Chenglin1，Yan Hongxiang2
（1.Department of Civil Engineering，University of Science and Technology Beijing，Beijing 100083，P.R.China；2.China Highway Engineering Consulting Corporation，Beijing 100097，P.R.China）

Based on the static liquefaction of granular materials，an elastoplastic static liquefaction constitutive model was proposed in the framework of critical state soil mechanics.The yielding surface and hardening rules were obtained by integrating the contact force of the model proposed by Chang and a nonassociate flow rule was adopted as well.The model has taken the state dependent dialatency law and the effect of initial density to the stress-strain relationship into consideration.The parameters of the model are simple and have certain physical meanings.The predicted results obtained by the model have a good agreement with the undrained triaxial test of Toyoura sand and sand-silt mixture.

granular material；static liquefaction；constitutive model；numerical simulation

TU443

A

1674-4764（2014）03-0011-07

10.11835／j.issn.1674-4764.2014.03.003

2013-09-10

国家自然科学基金（51178044）；新世纪优秀人才资助项目（NCET-11-0579）；中央高校基本科研业务费（FRFTP-12-001B）；北京高校“青年英才计划”

刘洋（1979-），男，博士，副教授，主要从事岩土工程研究，（E-mail）imaginationly＠163.com。