基于达维坚科夫骨架曲线的软土非线性动力本构模型研究

张如林，楼梦麟

（1.同济大学 土木工程防灾国家重点实验室，上海 200092；2.中国石油大学（华东） 储运与建筑工程学院，山东 青岛 266555）

1 引 言

土的动力非线性本构关系是表征土体动力特性的最基本关系，是土动力学研究的核心问题之一，同时也是岩土工程领域研究的热点和难点。模型建立得合理与否，对于重要工程场地的非线性地震反应，以及对于土工建筑及地基的动力响应机制的研究都是至关重要的。室内试验和现场观测均表明，土体的动力性质非常复杂，在循环荷载作用下，岩土材料的动应力-应变关系主要表现出非线性、滞后性和变形累积三方面的特征。在求解岩土体的地震反应问题时，需要考虑岩土体在动力作用下的非线性特征，就必然要涉及到岩土体的动力本构模型的选择问题[1]。

土体动力非线性本构模型研究主要分 3类方法：①工程中常用的等效线性化方法[2]，该方法形式直观简单，是试验结果的归纳，而且积累了大量工程经验，在目前土体动力反应分析、场地非线性地震反应分析中应用广泛，但它仍有不少缺点，如易产生“虚共振现象”，不能计算永久变形，在大应变、强震动时误差大，把土的塑性耗能转化为黏滞阻尼处理等等；②真非线性动力本构模型[3]，这类方法往往具有较多的模型参数，并且难以通过常规试验获得；③时域滞回非线性分析方法，该类方法包括Masing类、Iwan类以及多线性本构模型等，这也是本文研究所采用的方法。

试验研究表明，饱和软土的动剪切模量Gd与动剪应变γd之间的关系基本符合双曲线关系[4]。Hardin等[5]提出了预测Gd与γd的经验关系，即著名的Hardin-Drnevich模型，随后Martin等[6]认为，达维坚科夫模型可以更好地描述软土的动力变形特性。陈国兴等[7]对达维坚科夫骨架曲线进行了一定的修正，并与南京及其邻近地区土体试验结果进行对比，获得了达维坚科夫模型参数的拟合值。杨林德等[8]、刘齐建[9]曾对上海地区软土的动力特性进行了研究，通过对试验数据的拟合分析，发现对于上海典型饱和软土、上海软土的动力变形特性符合“应变软化”规律，认为达维坚科夫模型能够很好地描述Gd与γd、阻尼比D与γd之间的关系。

本文在前人研究成果基础上，对软土非线性动力本构关系进行了研究。基于Martin等[6]提出的达维坚科夫骨架曲线的土体动应力-应变关系，通过构造土体加载、再卸载的应力-应变关系滞回曲线，试图建立一种更为适合上海地区软土的非线性动力本构模型，从而为上海地区软土场地的动力计算反应分析提供适用的动力本构模型。

2 达维坚科夫骨架曲线和动应力-应变滞回曲线表达式

岩土体一维动应力-应变关系的研究是岩土材料动力本构关系研究中最基本的内容。Masing[10]最早建立了等幅循环荷载作用下岩土体的一维动应力-应变关系，按下列规则来构造动本构关系：①土的动应力-应变关系的骨架曲线为双曲线；②在初始加载过程中，动应力-应变关系遵循骨架曲线；③在初始反向卸载时土的动剪切模量与最大剪切模量相等，加、卸载应力-应变关系曲线与骨架曲线成2倍关系。上述规则被称为Masing法则，本节首先从一维应力-应变关系骨架曲线开始推导。

2.1 达维坚科夫骨架曲线理论公式

Hardin等[5]提出，用一种双曲线形式来描述土体的应力-应变关系，其动剪切模量比的表达式为

式中：G为剪切模量；Gmax为初始最大剪切模量；γ为剪应变幅值；γref=τfGmax，为参考剪应变，一般可取动剪切模量比为0.5所对应的剪应变幅值，如图1所示。

图1 双曲线应力-应变关系和参考剪应变示意Fig.1 Hyperbolic stress-strain curve and reference shear strain

Martin等[6]采用达维坚科夫骨架曲线来描述上述关系，将式（2）中的H(γ)改写成如下形式：

式中：A和B为与土性有关的拟合参数。需要指出的是，这里的γ0不再是具有明确物理意义的参考剪应变幅值，而仅仅是个拟合参数而已。显然，当A=1.0，B=0.5，γ0=γref时，达维坚科夫模型描述的应力-应变关系曲线就退化为Masing双曲线模型。

达维坚科夫模型应力-应变关系的骨架曲线可表示为

把式（3）代入到式（4），可得

2.2 骨架曲线增量剪切模量表达式

假定土体为理想黏弹性体，可从一般的模量衰减曲线得到归一化的剪应力为

式中：Ms为模量衰减系数，是一个归一化的割线模量。

通过对式（6）求γ的导数，可得到归一化的切线模量Mt，本文将其称为模量乘子，计算表达式如下：

增量剪切模量G就可以表示为

再对式（5）求导，就可得到达维坚科夫骨架曲线中的增量剪切模量表达式：

2.3 滞回曲线增量剪切模量表达式

当应力-应变关系的骨架曲线为达维坚科夫模型时，可采用Masing法则来构造相应的滞回曲线。根据Masing法则得到的滞回曲线的计算公式为

式中：τr为反转点的应力；γr为反转点的应变，如图2所示。

图2 骨架曲线和滞回圈的构造Fig.2 Skeleton curve and stress-strain hysteresis curves

联合式（10）和式（5）得到滞回曲线的计算表达式为

式中：R=(γ-γr)(2γ0)。

对式（11）求导，得出滞回曲线的增量剪切模量计算表达式为

至此，就得到了完整的达维坚科夫骨架曲线和卸载、再加载的滞回曲线的表达形式。

3 程序编制与验证

3.1 FLAC3D中本构模型的开发环境

本文在 FLAC3D提供的二次开发平台上进行前面的本构模型程序的编制。FLAC3D采用面向对象的标准C++语言编写而成，其所有的本构模型均以动态链接库文件（.DLL）的形式提供给用户，在计算过程中主程序会自动调用用户指定的本构模型的动态链接库文件[11]。FLAC3D自定义本构模型的主要功能是根据给出的应变增量得到新的应力。模型文件的编写主要包括 5部分内容：①基类（class constitutive model）的描述；②成员函数的描述；③模型的注册；④模型与FLAC3D之间的信息交换；⑤模型状态指示器的描述。由于 FLAC3D自带的本构模型和用户自己编写的本构模型继承的都是同一个基类，因此，用户自定义的本构模型和软件自带的本构模型的执行效率处于同一水平[11]。

对于本文来讲，FLAC3D的本构模型开发工作主要是修改头文件（.H文件）和程序文件（.CPP文件）。本文在线弹性本构模型（Userelas.H和Userelas.CPP两个文件）的开放源代码的基础上进行修改。在线弹性本构模型中，加载和卸载、再加载过程中的模量值一直保持恒定不变，而本文模型在每次计算过程中，都要根据当前的应变状态修改其模量值，来体现加载和卸载剪切模量的变化。

3.2 应变反转与实现思路

多数材料在一维条件下的力学性质与三维条件下的力学性质存在很大差别，尤其是对于岩土这样的颗粒破碎材料。在进行三维土层场地地震反应以及土与结构相互作用分析时，需要将前面的应力-应变关系从一维应变空间扩展到三维应变空间中。

在构建上述本构模型时，需要考虑应变反转的问题。本文采用如下处理方式，对于三维分析问题，至少存在6个应变率张量分量，分别定义如下[12]：

式中：Δeij为应变增量张量。

在应变空间中寻找应变轨迹的极值点，用上标0定义前一个点，用上标00定义前面第2个点，应变空间中前一个单位矢量的计算公式为

式中：下标i为1～6个分量。

新的应变增量 (εi-)在前一个单位矢量上的投影，或者用当前的偏移向量乘以就可以得到应变幅值：

将式（15）中的γ代入到式（7）中就可以得到模量乘子Mt，再代入到式（8）就可得到此时的剪切模量。

3.3 程序编制

Rosenblueth[13]、Newmark[14]等对 Masing 法则又做了2条补充规则：①若卸载应力-应变曲线或加载应力-应变曲线与骨架曲线相交，则偏离原曲线而遵循骨架曲线，称为“上骨架曲线”规则；②若当前曲线与先前同方向的后继应力-应变曲线相交，则当前曲线偏离原曲线，而遵循先前的后继应力-应变关系曲线移动，称“上大圈”规则。Masing法则和上述2条补充规则通常称为“扩展的Masing法则”，如图3所示。

图3 “上大圈”准则示意图Fig.3 Sketch of up-large-cycle curve rule

扩展的 Masing法则难以用简单的数学表达式描述，确定加、卸载的过程中应力-应变曲线的分支走向时只能采用数值判断来实现，在程序编制时需要特别对拐点加以判断。本文采用一种“栈”的数据结构形式，该数据结构采用“后进先出”的顺序来对其中的数据进行操作[15]。在执行过程中，每遇到一个拐点，便将其状态信息存入栈中（称为入栈），一旦当前的应变水平达到或者超过栈中记录的上次拐点的应变水平，就弹出记录的该拐点的状态信息（称为弹栈），进而更新剪切模量乘子，这样就可以简便地实现复杂加载路径下的应变转折。

FLAC3D自带一个Run()函数，该函数的主要作用是根据应变增量计算应力张量。本文在编制程序时，定义一个MyMt()函数来返回每个计算时刻的剪切模量值。在该函数中，首先提取Run()函数中的应变值，判断此时是否为拐点，以及是否超过栈中记录的上个拐点应变等操作，计算此时的模量乘子，进而更新模量值G并返回到Run()函数中，计算当前的应力值，如此重复这样的步骤直到计算结束。整个程序计算流程图如图4所示。

图4 程序计算流程图Fig.4 Flow chart of program

3.4 程序验证

在程序验证时，常采用一个单元进行计算，因为一个单元本身的应力和应变状态简单。为了验证程序编制的可靠性，检验程序能否实现应变拐点判断以及遵循广义 Masing法则等功能，以一个单元（1 m×1 m×1 m）为例，通过施加剪应变的方式来研究复杂加载路径下的应力-应变发展情况。剪应变加载路径如图5所示，计算得到的应力-应变关系发展曲线如图6所示。

可以看出，在复杂加载路径下，应变一旦遇到原来的拐点，仍会沿着原来的路径继续发展，当遇到骨架曲线，同样也会沿着骨架曲线继续前行，这说明本文编制的程序是合理可靠的。

图5 复杂加载路径曲线Fig.5 Curve of complex loading path

图6 计算的应力-应变关系曲线Fig.6 Calculating stress-stress curve by program

4 土体试验结果对比

4.1 相关文献土体模型参数

杨林德[8]、刘齐建[9]等通过对上海市区典型软土的动力试验，测定了试样的动剪切模量Gd、阻尼比D随动剪应变γd变化的关系曲线。本文采用文献[9]中粉质黏土典型土体试样的动剪切模量比GGmax-γd和D-γd的实测曲线，所选土样取自上海市轨道交通6号线德平路车站和明珠二期工程临平路车站，共有14个土样，其中德平路12个，临平路2个。根据文献[9]中试验数据点拟合的达维坚科夫模型参数，以及本文通过插值，取GGmax=0.5所对应的剪应变幅值作为 Hardin-Drnevich模型中的参考剪应变γr。

将本文结果与文献[16]中的土体试验结果进行对比，通过最小二乘法拟合出文献[16]中土体试样的达维坚科夫模型参数，以及数据插值得到Hardin-Drnevich模型的参考剪应变γr。文献[9]、[16]中两种土体试样的模型计算参数分别列于表1。

表1 两种本构模型计算参数Table1 Calculation parameters of two constitutive models

4.2 结果对比分析

根据表1所列的两种本构模型的计算参数，本文以一个单元（1 m×1 m×1 m）为例来进行计算分析。通过施加若干个水平的剪应变，得到不同剪应变幅值下动剪切模量比和阻尼比随剪应变幅值的变化曲线，并与文献[9]、[16]中的土样试验数据结果进行对比分析。对比结果如图7～10所示。

图7 计算G/Gmax-γd曲线与文献[9]试验数据对比Fig.7 Comparisons of simulated G/Gmax-γdcurves and test data of Ref. [9]

图8 计算D-γd曲线与文献[9]试验数据对比Fig.8 Comparisons of simulated D-γdcurves and test data of Ref. [9]

图9 计算G/Gmax-γ 曲线与文献[16]试验曲线对比Fig.9 Comparisons of simulated G/Gmax-γ curves and test data of Ref. [16]

图10 计算D-γ 曲线与文献[16]试验曲线对比Fig.10 Comparisons of simulated D-γ curves and test data of Ref. [16]

从图7～10中本文计算结果与文献[9]、[16]中的试验数据结果对比来看，在整个计算剪应变范围内，两种本构模型计算的动剪切模量比随剪应变幅值的变化曲线与试验曲线都拟合较好。相比Hardin-Drnevich模型，采用达维坚科夫骨架曲线构造的滞回曲线模型可以更好地拟合试验数据结果。在应变值小于 0.05%的范围内，达维坚科夫本构模型拟合的阻尼比曲线较好；在应变值大于 0.05%时，虽然Hardin-Drnevich模型和达维坚科夫模型拟合阻尼比曲线都有一定的差距，但达维坚科夫模型的拟合效果仍然比Hardin-Drnevich模型要好得多。

研究表明[17]，实际土层场地地震反应分析中，土体的动剪应变幅值一般不会超过0.1%，因此，在大部分的动剪应变范围内，采用达维坚科夫骨架曲线构造的非线性动力本构模型是可行的。从这个意义上来讲，本文所构建的基于达维坚科夫骨架曲线的非线性动力本构模型比传统的 Hardin-Drnevich模型要优越很多，也更适合上海地区软土场地的动力计算分析。

5 结 论

（1）上海地区软土的动力变形特性符合“应变软化”规律，可以用达维坚科夫模型描述。

（2）以达维坚科夫骨架曲线为基础，采用Masing法则，构造了土体加载、再卸载的动应力-应变关系滞回曲线，推导了骨架曲线和卸载、再加载滞回曲线的增量剪切剪切模量。

（3）将土体动应力-应变关系曲线从一维推广到三维应变空间，基于FLAC3D提供的二次开发平台，编制了基于达维坚科夫骨架曲线和符合广义Masing法则的土体非线性动力本构模型计算程序，通过复杂加载路径验证了编制程序的正确性。

（4）研究结果表明，相比工程上广泛应用的Hardin-Drnevich模型，基于达维坚科夫骨架曲线构造的本构模型所得到的软土 G/Gmax-γd和D-γd曲线更符合试验结果，验证了本文所建立的本构模型的合理性，可用于上海地区软土场地的动力计算反应分析。

（5）根据本文研究，在应变值超过 0.05%时，无论采用Hardin-Drnevich模型，还是达维坚科夫模型都会过高地估算阻尼比，由此而引起的计算误差还需做深入的研究。另外，软土动力特性十分复杂，影响软土土体非线性的因素也很多，特别是孔隙水压力效应、动态残余应变及残余应力特性、累积损伤及疲劳效应、震动软化等，有必要在今后的研究中加以适当考虑。

[1]庄海洋,陈国兴,梁艳仙,等. 土体动非线性黏弹性模型及其ABAQUS软件的实现[J]. 岩土力学,2007,28(3):1267－1272.ZHUANG Hai-yang,CHEN Guo-xing,LIANG Yan-xian,et al. A developed dynamic viscoelastic constitutive relations of soil and implemented by ABAQUS software[J]. Rock and Soil Mechanics,2007,28(3):1267－1272.

[2]齐文浩,薄景山. 土层地震反应等效线性化方法综述[J].世界地震工程,2007,23(4): 221－226.QI Wen-hao,BO Jing-shan. Summarization on equivalent linear method of seismic responses for soil layers[J].World Earthquake Engineering,2007,23(4): 221－226.

[3]WANG Z L,DAFALIAS Y F,SHEN C K. Bounding surface hypoplasticity model for sand[J]. Journal of Engineering Mechanics,ASCE,1990,116(5): 983－1001.

[4]HARDIN B O,RICHART F E. Elastic wave velocities in granular soils[J]. Journal of the Soil Mechanics and Foundations Division,ASCE,1963,89 (SM1): 33－65.

[5]HARDIN B O,DRNEVICH V P. Shear modulus and damping in soil: design equations and curves[J]. Journal of the Soil Mechanics and Foundation Engineering Division,ASCE,1972,98(7): 667－692.

[6]MARTIN P P,SEED H B. One dimensional dynamic ground response analysis[J]. Journal of Geotechnical Engineering,ASCE,1982,108(7): 935－954.

[7]陈国兴,庄海洋. 基于 Davidenkov骨架曲线的土体动力本构关系及其参数研究[J]. 岩土工程学报,2005,27(8): 860－864.CHEN Guo-xing,ZHUANG Hai-yang. Developed nonlinear dynamic constitutive relations of soils based onDavidenkov skeleton curve[J]. Chinese Journal of Geotechnical Engineering,2005,27(8): 860－864.

[8]杨林德,陆忠良,白延辉,等. 上海地铁车站抗震设计方法研究[R]. 上海: 同济大学防灾救灾研究所,2002.

[9]刘齐建. 软土地铁建筑结构抗震设计计算理论的研究[博士学位论文D]. 上海: 同济大学,2005.

[10]MASING G. Eigenspannungen und Verfestingung beim Messing[C]//Proceedings of the 2nd International Congress of Applied Mechanics. Zurich: [s.n.],1926: 332－335.

[11]陈育民,刘汉龙. 邓肯-张本构模型在 FLAC3D中的开发与实现[J]. 岩土力学,2007,28(10): 2123－2126.CHEN Yu-min,LIU Han-long. Development and implementation of Duncan-Chang constitutive model in FLAC3D[J]. Rock and Soil Mechanics,2007,28(10):2123－2126.

[12]CUNDALL P A. A simple hysteretic damping formulation for dynamic continuum simulations[C]//Proceedings of the 4th International FLAC Symposium on Numerical Modeling in Geomechanics. Minneapolis: Itasca Consulting Group,2006.

[13]ROSENBLUETH E,HERREAR I. On a kind of hysteretic damping[J]. Journal of the Engineering Mechanics Division,ASCE,1964,90(4): 37－47.

[14]NEWMARK N M,ROSENBLUTH E A. Fundamentals of earthquake engineering[M]. Englewood Cliffs: Prentice Hall Inc.,1971: 163－192.

[15]朱战立. 数据结构(C++语言描述)[M]. 北京: 高等教育出版社,2004.

[16]IDRISS I M,SUN J I. User's Manual for Shake91 Center for Geotechnical Modeling[D]. Davis: Department of Civil & Environmental Engineering,University of California,1992.

[17]CAMILO PHILLIPS,YOUSSEF M A. HASHASH.Damping formulation for nonlinear 1D site response analyses[J]. Soil Dynamics and Earthquake Engineering,2009,29(7): 1143－1158.