砂土液化大变形本构模型的三维化及其数值实现①

王 睿，张建民，王 刚

（1．清华大学水沙科学与水利水电工程国家重点实验室，北京 100084；2．清华大学土木水利学院，北京 100084；3．二滩水电开发有限责任公司，四川 成都 610021）

0 引言

砂土的液化后大变形作为地震液化诱发灾害的重要成因，自几次强震以来［1－2］一直是岩土工程抗震方面的一个研究热点。针对砂土在循环加载条件下的动力响应，很多学者进行了卓有成效的研究，并提出了一系列 包 括 边 界 面 模 型［3－5］、多 面 模 型［6－7］和 广义塑性模型［8］等砂土循环本构模型。但是这些模型中，鲜有能够基于砂土液化后大变形的机理有效模拟砂土液化后大变形的产生与累积。

基于试验观察，Shamoto和张建民［9］以及张建民和王刚［10－11］提出了饱和砂土液化后大变形的机理解释。将体变分解为有效球应力变化引起的体变、可逆性剪切体变和不可逆性剪切体变，揭示了不排水循环加载条件下砂土液化后剪应变的发展和累积与三个体积应变分量之间的内在联系。根据体积相容性条件，定量的描述了砂土液化后在零有效应力条件下所产生的“似流体”剪应变 （fluid－like shear strain）γ0。

以此机理为基础，王刚和张建民［12］在边界面塑性理论［13］框架下，建立了一个可以描述饱和砂土液化后大变形的弹塑性循环本构模型。该模型根据试验规律定义了区别于其他已有模型的可逆和不可逆剪胀速率，从而实现了对砂土循环剪切体变的正确描述。通过体积相容性条件，有效地计算砂土液化后剪切变形的产生与累积，为定量描述砂土液化后大变形提供了一个科学的建模方法。由于三维条件下可逆性剪胀面上的应力映射规则并未给定，以及根据边界面上映射规则确定的流动方向数值上求解成本过大，该模型可以在二维应力空间中应用，但无法用于真正的三维计算。

基于该模型，本文提出了与其相适应的三维应力空间中边界面和剪胀面上的应力映射规则，并根据数值计算能力需求，调整了塑性流动方向，对砂土液化大变形本构模型进行了三维化。利用高效的数值算法，在OpenSees［14］开源有限元平台上对模型进行三维的数值实现；通过对试验的模拟和三维地基动力反应分析，验证模型和数值算法的有效性。

1 模型的三维形式

由于王刚和张建民［11－12］对本模型大部分要素做了详细的描述，本文重点讨论在三维应力空间中模型的相关映射规则，其他方面仅简要介绍。

本文各处张量均采用加黑斜体，标量采用斜体表示。除特别声明变量外，变量大多采用土力学常规符号表达。

1．1 基本应力应变关系

模型采用弹塑性应力应变关系

1．2 弹性模量

弹性模量采用Richart［15］建议的方式定义

pa为标准大气压。

1．3 破坏面、边界面和可逆性剪胀面

由于本模型为边界面塑性理论［13］〛框架下的本构模型，在确定各种映射规则前，首先需要确定本模型所使用的破坏面、边界面和可逆性剪胀面，分别定义如下：

式中Mf，c为三轴压缩条件下的破坏应力比；ηm，c为历史最大三轴压缩应力比；Md，c为三轴压缩条件下可逆性剪胀的产生与释放的分界；θ为应力罗德角；而形状函数g（θ）则为根据张建民［16］的简化函数（图1）。

1．4 加载与反向加载判断

模型通过塑性加载指数符号判定加载与反向加载：

1．5 边界面上的映射规则与塑性模量

在边界面塑性理论框架下［13］，塑性模量受到当前应力点到其在边界面上投影的距离控制，因此需要定义恰当的映射规则。

对于当前应力在边界面上的投影¯r，借鉴王志良等［3］提出的相应映射规则，定义为上一次反向加载点（应力反转点）ain到当前应力点r连线的延长线与边界面的交点（图1）。因此可以将边界面上投影点表示为

其中β可以通过将式（9）带入式（6）求解。

图1 模型映射规则Fig．1 Mapping rules for the model．

假定偏应力空间中的加载方向n与偏应变增量方向m一致，并定义为沿¯r方向的单位张量：

实际上，严格根据弹塑性理论n应当为边界面在应力投影点处的单位外法向，但是由于在模型数值化时，对于边界面函数求解外法向在数值上十分困难且计算成本很高，因此在模型的三维化中采用这一调整后的加载与流动方向。

在定义了模型在边界面上的应力投影规则后，塑性模量的确定可以通过下式得到：

1．6 可逆剪胀面上的映射规则和可逆剪胀速率

由于王刚和张建民［12］并未定义三维空间中可逆剪胀面上的应力投影规则，若采用模型二维形式中可逆性剪胀生成与释放的判别标准和剪胀速率，将在三维空间中出现误差。以图1中的应力状态为例，按照模型二维形式中的判断规则，由于｜η｜≥Md，cg（θ）且｜η｜＞0，因此判断为发生可逆性剪胀，可逆性剪胀速率Dre＜0。而实际上对于三维应力空间，可逆性剪胀／减缩的判断应当由当前应力点r在应力反转点ain到r的方向上与可逆性剪胀面的相对位置确定。由此，定义当前应力点r在可逆性剪胀面上的投影rd为（图1）：

在该映射规则下，可以定义三维应力空间中可逆性剪胀的产生与释放判定标准为

式中，Dre，gen为可逆性剪胀产生速率；Dre，rel为释放速率。

相应的，三维条件下可逆性剪胀产生速率为

可逆性剪胀释放速率仍为［12］

dre，2为可逆性剪胀释放参数。

根据新的映射规则和可逆性剪胀的产生与释放判定标准，图1所示情况发生可逆性剪缩，可逆性剪胀速率Dre＝Dre，rel＞0，这是由于当前应力点处于其所对应的边界面上投影点“以内”。

1．7 不可逆剪胀速率

不可逆剪胀速率仍参照王刚和张建民［12］：

dir，α和γd，r为不可逆性剪胀参数；γmono为单向剪应变长度。

1．8 液化后大变形

通过将式（4）带入式（1）积分，可以得球应力发展至零时的εvc阈值，即εvc，0：

当εvc大于εvc，0时，有效球应力保持为零，土体进入液化状态，式（1）中的体变公式不再有效。εvc转而由体积相容性确定：

此时各剪胀关系依然有效，因此只有当足够的剪应变ep发生时（即液化后大变形），才能产生充分的εvd，re，使得土体离开液化状态［10－11］。

2 模型数值实现

通过高效数值算法，在OpenSees［14］中进行了三维化模型的数值实现。

2．1 零有效应力处理方法

关于零有效应力状态下模型的数值处理，采用王刚和张建民［12］提出的方法，即设定一最小有效球应力pmin，使得

此时，εvc，0可表示为

2．2 应力积分算法

采用分子步的Cutting Plane半显式积分算法进行应力积分。对于半显式积分算法，需要控制每一步应力积分时的应变步长，否则可能出现不收敛的情况。为了增加算法稳定性，将有限元得到的应变增量分成子步后再进行应力积分计算。

计算子步长时，首先对下一步应力进行弹性预测，当剪应变增量或剪应力比增量超过容许值，则将步长拆分为n个子步：

对于每一个子步采用Cutting Plane算法进行应力积分。该算法首先由Simo和Ortiz［17］提出，主要思路是在弹性预测应力增量的基础上，通过对一致性条件φ＝0中函数φ的一阶展开逼近一致性条件进行塑性修正，最后得到实际应力增量，步骤如下：

（1）初始化局部迭代次数、塑性应变和塑性加载指数：

（2）弹性预测：

（3）检查一致性条件，判断加载／反向加载：

当φ（k）＞0，发生塑性加载，进行第4步；当φ（k）＜0，应力反转，更新应力反转点（ain）n＋1＝rn，进行第6步。

（4）通过计算塑性加载指数增量进行塑性修正：

更新塑性加载指数和应力应变状态：

（5）检查一致性条件残余值：

若｜φ（k＋1）｜＞tolerance，则未收敛，k＝k＋1，进行步骤4；否则进行第6步。

（6）更新应力和内变量：

对于半显式的Cutting Plane积分算法，其四阶刚度张量如下：

其中De为弹性刚度张量。

2．3 边界面上应力投影点求解算法

在应力积分的第3和第4步中，需要求解应力在边界面上的投影，即求解式9中的β。为此采用Dowel和Jarratt［18］提出的Pegasus求根算法，该算法能够保证快速无条件收敛，具体步骤如下：

（1）初始化试探值β0＝0和β1＝1。

（2）计算

（4）计算

若｜fb（β）｜＜tolerance，计算收敛；否则进行第5步。

通过上述算法，目前已经在 OpenSees 2．4．0版本中加入了三维化的模型。为了增加初始静应力场计算时模型的数值稳定性，添加了弹性材料阶段，具体使用方法可参见OpenSees网站中该模型说明文档（http：／／opensees．berkeley．edu）。

3 三维化模型及算法的应用

利用OpenSees中已有的用基于饱和土动力固结理论的完全耦合单元，使用本模型对饱和砂土不排水循环扭剪试验进行模拟，并通过一个三维倾斜地基动力反应分析验证模型及相应数值算法对于计算液化后大变形和处理三维问题的有效性。

3．1 内华达砂不排水循环扭剪试验模拟

对Kutter等［19］用内华达砂进行的不排水空心圆柱扭剪试验进行模拟。内华达砂的基本物性指标为：Gs＝2．67；d50＝0．18mm；emax＝0．887；emin＝0．551。试验中所用砂土相对密度71%，采用等剪应力幅值为56kPa进行循环加载，围压198kPa。由于该试验目的之一是用来验证王志良等［3］的模型和其他边界面模型的有效性，因此弹性和塑性模量以及破坏参数均可由Kutter等［19］的试验报告中直接获得；剪胀相关参数则采用相应的参数确定方法估计。具体参见表1。

图2 不同降雨量时水平位移等值线Fig．2 Horizontal displacement field at different rainfall level．

表1 模型参数取值Table 1 Parameter values in the model

图2给出了对该试验的模拟结果，可以看到计算结果与试验结果的应力路径和应力应变关系吻合较好，尤其是对初始液化后应力路径和零有效应力时产生的液化后大变形的模拟非常有效。

3．2 三维倾斜地基动力反应分析

利用OpenSees中Brick UP［20］三维完全耦合单元，对一个10m倾斜地基进行了三维动力反应分析。如图3所示，地基向y方向倾斜2°，地震动输入沿x方向和z方向双向输入。模型网格由10个边长1m的立方体单元组成，底面节点约束y向位移；为了模拟无限地基，将同层结点的三向自由度绑定；模型底面和侧面为不排水面，地基表面排水，孔压保持为零。土体采用VELACS离心模型试验［20］中使用的相对密度为45%的内华达砂的参数。输入地震波采用VELACS离心模型试验［21］模型1所使用的波形，竖向振幅设定为水平向的0．2倍。

图3 模型示意图Fig．3 Schematic of model set up．

图4给出了地震全过程不同时刻土体的变形和孔隙水压力云图。可以看出，在最初的2～3s内土体尚未达到液化，地基变形较小。当浅层土体孔压累积至液化，地基开始向y方向变形，直至地震结束后在1m深处达到0．13m。图5给出了1m深度处孔压、x和y方向位移时程。超静孔压比在4s左右初次达到1，即初始液化。x方向位移随着地震波而波动，但地震结束后残余值几乎为零；y方向位移则随着超静孔压比接近1迅速累积。

对这一倾斜地基的三维动力反应分析表明，本模型和相应的数值算法能够有效地计算三维条件下砂土的动力响应和液化变形。

图4 地震过程中土体变形（放大10倍）和孔压云图Fig．4 Soil deformation（amplified to 10times）and pore pressure contour during earthquake．

图5 1m深度处孔压及x和y方向位移时程Fig．5 Pore pressure and xand ydirection displacement history at 1mdepth．

4 结论

本文对基于符合砂土液化后大变形机理的弹塑性本构框架，建立了三维应力空间中砂土液化大变形本构模型，实现了其数值化，并用于试验模拟和三维动力反应分析。

发展了在三维应力空间中适合砂土液化后大变形本构模型的边界面和剪胀面上应力映射规则，并相应的定义了新的可逆性剪胀产生率。根据数值计算能力需求，调整了塑性流动方向和加载方向。实现了砂土液化后大变形本构模型的三维化。

通过采用分子步的Cutting Plane数值积分算法实现模型应力积分，并利用高效无条件收敛的Pegasus求根算法求解应力投影，在OpenSees开源有限元平台上实现了三维模型的数值化，对所有研究人员开放使用。

利用该模型，对内华达砂不排水循环扭剪试验进行模拟，验证了模型的模拟能力，尤其是在计算砂土液化后大变形方面的优势。通过对一个三维倾斜地基的动力反应分析，显示了模型及相应数值算法在处理三维问题上的有效性。

［1］ Seed H B．Soil liquefaction and cyclic mobility evaluation for level ground during earthquakes［J］．Journal of the Geotechnical Engineering Division，1979：105（2）：201－255．

［2］ Hamada M．Large Ground Deformations and their Effects on Lifelines：1964Niigata Earthquake．Case Studies of Liquefaction and Lifelines Performance during Past Earthquake［R］．Buffalo：National Centre for Earthquake Engineering Research，1992：NCEER－92－0001．

［3］ Wang Z L，Dafalias Y F，Shen C K．Bounding surface hypoplasticity model for sand［J］．Journal of Engineering Mechanics，1990，116（5）：983－1001．

［4］ Papadimitriou A G，Bouckovalas G D，Dafalias Y F．Plasticity model for sand under small and large cyclic strains［J］．Journal of Geotechnical and Geoenvironmental Engineering，2001，127（11）：973－983．

［5］ Dafalias Y F，Manzari M T．Simple plasticity sand model accounting for fabric change effects［J］．Journal of Engineering Mechanics，2004，130（6）：622－634．

［6］ Mroz Z，Norris V A，ZienkiewtczОC．An anisotropic hardening model for soils and its application to cyclic loading［J］．International Journal for Numerical and Analytical Methods in Geomechanics，1978，3（2）：203－221．

［7］ Yang Z，Elgamal A，Parra E．Computational model for cyclic mobility and associated shear deformation［J］．Journal of Geotechnical and Geoenvironmental Engineering，2003，129（12）：1119－1127．

［8］ Pastor M，Zienkiewicz O C，Chan A．Generalized plasticity and the modelling of soil behavior［J］．International Journal for Numerical and Analytical Methods in Geomechanics，1990，14（3）：151－190．

［9］ Shamoto Y，Zhang J M．Mechanism of large post－liquefaction deformation in saturated sands［J］．Soils and Foundations，1997，2（37）：71－80．

［10］ 张建民，王刚．砂土液化后大变形的机理［J］．岩土工程学报，2006，28（7）：835－840．

Zhang J M，Wang G．Mechanism of large post－liquefaction deformation in saturated sand［J］．Chinese Journal of Geotechnical Engineering，2006，28（7）：835－840．

［11］ Zhang J M，Wang G．Large post－liquefaction deformation of sand，part I：physical mechanism，constitutive description and numerical algorithm［J］．Acta Geotechnica，2012，7（2）：69－113．

［12］ 王刚，张建民．砂土液化大变形的弹塑性循环本构模型［J］．岩土工程学报，2007，29（1）：51－59．

Wang G，Zhang J M．A cyclic elasto－plastic constitutive model for evaluating large liquefaction－induced deformation of sand［J］．Chinese Journal of Geotechnical Engineering，2007，29（1）：51－59．

［13］ Dafalias Y F，Popov E P．A model of nonlinearly hardening materials for complex loading［J］．Acta Mechanica，1975，21（3）：173－192．

［14］ McKenna F，Fenves G L．OpenSees Manual［EB／OL］．PEER Center，http：／／OpenSees．berkeley．edu．2001．

［15］ Richart F E，JR，HALL J R，Woods R D．Vibrations of soils and foundations［M］．Englewood Cliffs，N．J．：Prentice－Hall．Inc．，1970．

［16］ Zhang J M．Cyclic critical stress state theory of sand with its application to geotechnical problems［R］．Tokyo：Research Report of Tokyo Institute of Technology，1997：1－269．

［17］ Simo J C，Ortiz M．A unified approach to finite deformation elastoplastic analysis based on the use of hyperelastic constitutive equations［J］．Computer Methods in Applied Mechanics and Engineering，1985，49（2）：221－245．

［18］ Dowell M，Jarratt P．The"Pegasus"method for computing the root of an equation［J］．Bit Numerical Mathematics，1972，12（4）：503－508．

［19］ Kutter B L，Chen Y，Shen C K．Triaxial and torsional shear test results for sand［R］．CR 94．003－SHR，Port Hueneme，California：Naval Facilities Engineering Service Center，Contact Report，1994．

［20］ Yang Z，Lu J，Elgamal A．OpenSees soil models and solidfluid fully coupled elements user manual［M］．San Diego，California：University of California，2008．

［21］ Taboada V M，Dobry R．Experimental results of Model No．1at RPI［A］∥Proceedings of International Conference on Verification of Numerical Procedures for the Analysis of Soil Liquefaction Problems［C］．［S．l．］：［s．n．］，1993，3－17．