基于颗粒流程序的广义Kelvin 模型及其应用

金爱兵，王 凯，张秀凤，孟新秋，杨振伟

（1.北京科技大学 土木与环境工程学院，北京 100083；2.温州市城建设计院，浙江 温州 325000）

1 引 言

土体作为工程材料的一种，在外荷载作用下的力学响应不能完全由经典力学描述。特别是软土，由于其内部结构特殊，其应力-应变关系呈现出一定的非线性流变特征。陈宗基[1]通过将Maxwell 体与Hook 弹性体并联，得到一种新模型来描述土体变形规律；孙钧[2]对上海地区3 种典型软土进行了大量的单轴剪切与三轴压缩流变试验，并对试验数据进行详细分析后，提出了上海软土蠕变特性的方程式。粉质黏土作为一种常见软土，也常常表现出一定的流变特性。王坤等[3]对深圳地区深层粉质黏土进行了低围压作用下的三轴流变试验，并得到了粉质黏土的三轴流变规律。而对于服役时间较长的露天矿山边坡，粉质黏土的流变效应对边坡工程的安全影响更为显著。

目前，对岩土流变特性的研究主要集中在室内试验[4－6]和流变模型开发[7－9]等方面。室内试验研究主要有单轴试验[6]和不同围压下的三轴试验[10－14]。熊良宵等[15]将修正后的Burgers 模型应用于蠕变试验的拟合分析；李娜等[16]用改进的西原模型对循环加、卸载试验结果进行拟合，所得效果良好。但室内试验研究往往存在操作复杂、花费时间长、试验费用高等缺点。

作为岩土工程稳定性研究的常用工具，数值模拟方法已经被许多科研工作者应用于岩土流变特性的研究。张志沛等[17]利用FLAC3D中的Burgers 模型对不同形状的试样进行了蠕变特性模拟分析；熊良宵等[18]在FLAC3D中开发出新的流变模型，并将其应用到岩土应变软化的相关工程计算中；刘珊珊等[19]在FLAC3D中二次开发出黏弹性广义Kelvin 模型，并用算例验证了所开发模型的正确性；王涛等[20]开发出适用于PFC2D的广义Kelvin 接触本构模型，并将其应用到了隧洞工程问题中。

本文在对PFC2D中自定义本构模型进行研究的基础上，利用Microsoft Visual Studio 平台，开发出广义Kelvin 模型的动态链接库(DLL)文件。通过算例验证程序编制的准确性，并在此基础上，将其应用到露天煤矿含粉质黏土边坡的失稳研究中。

2 广义Kelvin 模型本构关系

广义Kelvin 模型由一个Maxwell 体和任意个Kelvin 体串联组成，其模型如图1 所示。其中，Maxwell 体由一个弹簧和一个阻尼器串联而成，Kelvin 体由一个弹簧和一个阻尼器并联而成，Kn和Ks代表弹簧法向和切向刚度，Cn和Cs代表阻尼器法向和切向阻尼，下标数字1,2…,n为原件编号。

图1 广义Kelvin 模型Fig.1 Generalized Kelvin model

由图1 可知，广义Kelvin 模型的总位移u 等于一个Maxwell 体的位移 um1与若干个Kelvin 体的位移ukr（下标r为开尔文体的编号，r=2,3,…,n）之和，如式（1）所示。其中，um1由Maxwell 体中弹簧的位移 uk1和阻尼器的位移 uc1共同组成。

对广义Kelvin 模型的本构关系推导如下，下列各式中的“±”(或“∓”)分别对应着法向和切向的情况；Kr和Cr(r=1,2,…,n) 分别为广义Kelvin 模型中弹簧的刚度和阻尼器的阻尼；f为模型中颗粒间的接触力；Ar、Br、C和D 是分别与Kr、Cr和Δt 有关的中间变量。

对式（1）进行一阶求导，得

接触力f 可表示为

将式（3）、（5）、（6）代入式（2）得

利用有限差分格式的中心差分法近似对时间t进行差分，其中对t 取平均值得

对于Kelvin 体，则有

整理后，接触力ft+1可表示为

ft+1可由已知量ut+1、ut、和ft求出。

3 广义Kelvin 模型开发

PFC 程序中的用户自定义接触模型的编写包括基类描述、模型注册、成员函数描述、模型和PFC程序间的信息交换。软件还提供了自带本构模型的源代码，用户在进行本构模型开发时，可以利用和修改这些代码。下面介绍如何修改本构模型的头文件和源文件，从而创建用户自己的本构模型。

在头文件中进行新的本构模型派生类的声明，修改模型的编号、名称和版本，修改派生类的私有成员，包括模型的基本变量及主要中间变量。为避免冲突，新自定义模型的编号需不小于100。

在源文件中修改以下各个成员函数：

（1）修改构造函数Kelvin::Kelvin()，对其私有成员进行初始化。

（2）修改Kelvin::Name()函数，返回用户自定义模型的名称。

（3）根据派生类私有成员的名称字符串，依次修改Kelvin::PropNames()、Kelvin::ReturnProp()和Kelvin::AcceptProp()3 个函数。PFC 的PROPERTY命令可用来调用Kelvin::AcceptProp()函数，对模型的参数进行赋值。

（4）修改Kelvin::FDlaw()函数，此函数在计算离散元的力与位移关系时被调用。用户根据推导的模型中力与位移数学关系，编写C++语句，由应变增量计算得到应力增量，从而获得新的应力。本构模型将在每个循环中更新颗粒间接触的力与力矩，以及法向和切向刚度。

（5）修改USERCM1::SaveRestore()函数，方法同步骤（3），此函数用来保存与恢复实型、整型和逻辑型的内部数据，PFC 程序给出save和restore命令时调用此函数。

另外，针对本文自定义的广义Kelvin本构模型，程序中设置控制变量i来灵活控制Kelvin体的个数。编写Kelvin.h和Kelvin.cpp 文件后，编译动态链接库文件。至此，所开发的广义Kelvin 模型以动态链接库文件的形式存储起来，可供PFC2D程序调用。

4 广义Kelvin 模型验证

当控制变量i=0 时，即广义Kelvin 模型中只含有Maxwell 体，无Kelvin 体时，模型即为Maxwell模型；当控制变量i=1 时，即广义Kelvin 模型中只含有Maxwell 体和1 个Kelvin 体时，模型即为Burgers 模型。本文分别采用i=0和1 两种情况下的Maxwell 模型和Burgers 模型进行松弛试验的验证。

（1）Maxwell 模型理论值求解

当i=0，无Kelvin 体时，式（7）可表示为

（2）Burgers 模型理论值求解

当i=1 时，对于单位输入函数 u( t)=1，Burgers模型的理论值为

两种模型中Maxwell 体弹簧法向、切向刚度：K1n=K1s=100MPa；阻尼器法向、切向阻尼：C1n=C1s=100MPa·s；Kelvin体弹簧法向、切向刚度：K2n=K2s=100MPa，阻尼器法向、切向阻尼：C2n=C2s=100MPa·s。

通过上述计算，得到两种接触模型松弛过程的理论值后，采用PFC2D命令流和Fish 语言编写应力松弛程序，然后在PFC2D软件中调用该程序，分别加载i为0和1 时的广义Kelvin 本构模型动态链接库文件，模拟两个相互接触球体（两球重叠量为0.01 m）的应力松弛过程，得出球体间法向接触力随时间的变化关系，并将数值模拟解与理论值进行对比分析。其对比结果如图2 所示。

由图2 可知，两球接触力随时间的增加呈指数趋势减小，本构模型的数值模拟解与理论解所对应曲线相互吻合，从而验证了所开发模型的准确性。

图2 法向接触力模拟计算结果与理论结果Fig.2 Analytical and numerical results of normal contact forces

5 工程应用

5.1 边坡失稳规律

含软弱粉质黏土边坡广泛分布于山西、陕西等黄土高原地区，尤其是在山西地区，煤炭资源丰富，露天开采形成的含粉质黏土边坡普遍存在。在山西某露天煤矿，经长期实地调研，发现此类边坡的滑坡形式较为特殊，其失稳规律有以下几个特征：①边坡的滑坡区域含有粉质黏土层；②岩土接触带倾向坡内，倾角较小(6°～9°)，且有地下水渗出；③滑坡体前缘平缓，局部有底鼓隆起现象；④滑坡体中部较为平缓；⑤滑坡体后缘较陡，张裂缝严重，下落量小；⑥ 滑动面形式表现异常，为“圆弧+折线”的组合型，最终轮廓表现为帽檐形。

此类含粉质黏土边坡失稳破坏的一般过程与现场情况如图3 所示。

5.2 模型计算参数

岩土体的物理力学参数由岩石压缩试验、土工直剪试验及粉质黏土三轴流变试验获得，见表1。参照室内试验，进行数值模拟试验，粉质黏土部分调用本文自主开发的广义Kelvin 模型，其他岩土体均采用PFC2D软件默认线弹性接触模型；数值模型尺寸均为100 mm×50 mm 的标准试样；岩土体颗粒间均采用能够抗拉、抗剪、承受弯矩的平行黏结。

模拟试验过程中，对广义Kelvin 模型中的控制变量i 及模型的细观力学参数进行反复调试。最终结果表明，当控制变量i=2，即广义Kelvin 模型中含1 个Maxwell 体和2 个Kelvin 体时，粉质黏土的模拟试验与室内试验结果较吻合（见图4、5）。此时，在广义Kelvin 模型中，Maxwell 体弹簧法向、切向刚度：K1n=K1s=100 MPa；阻尼器法向、切向阻尼：C1n=C1s=300 MPa·s；两个Kelvin 体相同，弹簧法向、切向刚度：K2n=K2s=K3n=K3s=100 MPa；阻尼器法向、切向阻尼：C2n=C2s=C3n=C3s=300 MPa·s 。边坡模型中岩土体的其他细观力学参数见表2。

图3 边坡失稳破坏过程与现场Fig.3 Process and scene of slope failure

表1 岩土体力学参数Table 1 Mechanical parameters of rock and soils

图4 数值试验与直剪试验结果Fig.4 Results of numerical test and direct shear test

图5 数值试验与流变试验结果Fig.5 Results of numerical test and rheological test

5.3 边坡模型构建

根据参数调试过程中获取的广义Kelvin模型和细观力学参数，建立边坡开挖前的PFC2D数值计算模型，见图6。根据工程实际情况，边坡模型长×高尺寸为：217.51 m×82.66 m。模型内部共生成36 261 个岩土体颗粒，模型的外部由1 道水平墙体和2 道竖直墙体构成边界限制条件，墙体法向刚度取 kn=1 000 MN/m，切向刚度取 ks=1 000 MN/m，摩擦系数取1.0。模型中不同颜色颗粒分别代表不同地层，边坡自上而下依次为粉土、粉质黏土（天然）、粉质黏土（饱和）、接触带及泥岩。

表2 颗粒细观力学参数Table 2 Meso mechanical parameters of particles

图6 边坡数值模型Fig.6 Numerical model of slope

对于粉质黏土，考虑到吸水-失水循环对其强度的影响，根据现场实际情况，在边坡模型的不同地层位置和计算的不同时段，对粉质黏土分别采用天然强度与饱和强度。无降水时下层粉质黏土采用天然强度，降水时下层粉质黏土采用饱和强度。后续数值模拟中，边坡开挖后通过在计算的不用时步对该层粉质黏土赋予相应参数来实现这一过程，根据当地雨季的降水规律：模拟中每40 万时步采用天然强度，每10 万时步采用饱和强度，循环多次至边坡最终破坏形成，共计约240 万时步。

5.4 计算结果及分析

在PFC2D中编写命令流，实现对边坡数值模型的计算与开挖。对于边坡中的粉质黏土地层，分别采用线弹性接触模型和本文所开发的具有流变特性的广义Kelvin 模型两种情况进行模拟计算，两算例中吸水-失水循环模拟过程相同。最终将其计算结果进行对比分析。

计算中采用颗粒间平行黏结的拉裂或剪断来表示边坡岩土体结构的破坏[21]，为便于观察，将模型中黏结破坏的颗粒设置为深蓝色。模型中粉质黏土地层采用线弹性接触模型的计算结果如图7 所示，边坡最终破坏模式呈常见的圆弧形土体滑坡模式，但与现场边坡实际破坏情况不符。

图7 线弹性接触模型边坡破坏结果Fig.7 Slope failure predicted by a linear elastic contact model

模型中粉质黏土地层采用用户自定义广义Kelvin 模型的计算结果，如图8 所示，边坡最终破坏模式比较符合现场的实际情况。计算过程中，边坡的失稳破坏可分为5 个阶段。

第1阶段：边坡施工开挖，如图8(a)所示，各台阶开挖临空面小范围内岩土体受到扰动，深蓝色颗粒开始出现，并零星散布于边坡面的浅层。

第2阶段：边坡局部破坏，如图8(b)所示，边坡内部整体变化微小，因赋存地下水，粉质黏土层强度降低，其坡脚处开始有深蓝色颗粒聚集，局部位移增大，破坏初显。

第3阶段：破坏范围扩大，如图8(c)所示，饱和粉质黏土层边坡坡脚处深蓝色颗粒增多，孔隙率增加，破坏范围扩大，少量黏结已被破坏的颗粒在重力作用下开始滑移出坡面，边坡后缘开始呈现出较明显的圆弧滑面。

第4阶段：边坡出现失稳，如图8(d)所示，坡面周边大量深蓝色颗粒翻滚下滑，孔隙率大幅增加，边坡破坏区域轮廓成“折线+圆弧”的帽檐形。

第5阶段：帽檐形滑坡形成，如图8(e)所示，模型中的颗粒停止下滑，边坡模型已经破坏。

图8 广义Kelvin 模型边坡破坏过程Fig.8 Slope failure process derived from the generalized Kelvin model

通过上述计算，并结合现场滑坡情况，边坡失稳过程主要表现为：

（1）天然状态下，粉质黏土层强度较高，所在边坡具有足够的自稳能力，随着开采的进行，粉质黏土层中颗粒黏结遭到破坏，导致裂隙产生，计算过程中坡体共产生498 条裂隙。

（2）微观裂隙贯通产生宏观裂缝，受水影响，饱和粉质黏土在接触带处强度较低，颗粒黏结容易遭到破坏，所以该处产生裂隙最多，成为边坡破坏的主要区域。

（3）受开采影响，失水-吸水循环过程中，裂隙进一步扩张，松散体影响范围逐渐扩大，并在边坡后缘开始形成圆弧滑面。

（4）因松散体自身承载力不够，内部土体颗粒开始不断下滑，在粉质黏土流变作用下，一段时间后，粉质黏土不能承受上部土体的荷载，松散体呈折线形式向上蔓延，联通至边坡后缘圆弧滑面，以“折线+圆弧”的组合滑动面形式破坏，破坏区最终轮廓呈帽檐形。

含粉质黏土边坡失稳模拟结果表明，该类边坡失稳与普通土质边坡的失稳形态以及失稳过程并不一致，粉质黏土的流变性对该类边坡失稳具有重要影响；而水是在粉质黏土流变过程中导致边坡失稳的重要原因，因此，实际工程中可以采取压坡阻水的方式处治该类边坡，以保证边坡稳定。

6 结 论

（1）自定义广义Kelvin 本构模型包含1 个Maxwell 体和若干个Kelvin 体，实际应用中可以根据具体情况选择合适数量的Kelvin 体，以便更好地模拟岩土体流变特性。

（2）在PFC2D软件中调用基于C++语言开发的自定义广义Kelvin 本构模型动态链接库，进行算例验证，得到的数值模拟解与理论解吻合，模型具有较好的准确性。

（3）利用PFC 调用自定义广义Kelvin 本构模型对国内某含粉质黏土边坡模拟表明，含1个Maxwell体和2 个Kelvin 体的用户自定义广义Kelvin 模型可很好地模拟粉质黏土的流变特性，得到的边坡失稳破坏过程与现场基本一致。

（4）由于流变特性的影响，该类粉质黏土边坡破坏形式较为特殊，与常规土质边坡破坏呈圆弧滑动不一致，其破坏区最终轮廓成“折线+圆弧”的帽檐形。

[1]陈宗基.固结及时间效应的单维问题[J].土木工程学报,1958,5(1):1－10.CHEN Zong-ji.Unidimensional problems of consolidation and time effect[J].China Civil Engineering Journal,1958,5(1):1－10.

[2]孙钧.岩土材料流变及其工程应用[M].北京:中国建筑工业出版社,1999.

[3]王坤,蔡永昌,魏赟,等.深圳地区深层粉质黏土的流变特性试验研究[J].岩土工程界,2009,12(12):86－89.WANG Kun,CAI Yong-chang,WEI Yun,et al.Experimental study on the rheological properties of deep silty clay in Shenzhen[J].Geotechnical Engineering World,2009,12(12):86－89.

[4]ZHANG ZHILIANG,XU WEIYA,WANG WEI,et al.Triaxial creep tests of rock from the compressive zone of dam foundation in Xiangjiaba Hydropower Station[J].International Journal of Rock Mechanics and Mining Sciences,2012,50(1):133－139.

[5]胡华,郑晓栩.动载作用频率对海相沉积软土动态流变特性影响试验研究[J].岩土力学,2013,34(增刊1):9－13.HU Hua,ZHENG Xiao-xu.Experimental research on dynamic rheological characteristics of marine deposit soft soil under different frequencies of dynamic loading[J].Rock and Soil Mechanics,2013,34(Supp.1):9－13.

[6]LI YONGSHENG,XIA CAICHU.Time-dependent tests on intact rocks in uniaxial compression[J].International Journal of Rock Mechanics and Mining Sciences,2000,37(3):467－475.

[7]ZHANG HUABIN,WANG ZHIYIN,ZHENG YALI,et al.Study on triaxial creep experiment and constitutive relation of different rock salts[J].Safety Science,2012,50(4):801－805.

[8]LIU LANG,WANG GEMIN,CHEN JIANHONG,et al.Creep experiment and rheological model of deep saturated rock[J].Transactions of Nonferrous Metals Society of China,2013,23(2):478－483.

[9]SHAO J F,ZHU Q Z,SU K.Modeling of creep in rock materials in terms of material degradation[J].Computers and Geotechnics,2003,30(7):549－555.

[10]ZHANG YU,XU WEIYA,GU JINJIAN,et al.Triaxial creep tests of weak sandstone from fracture zone of high dam foundation[J].Journal of Central South University,2013,20(9):2528－2536.

[11]YANG SHENGQI,JIANG YUZHOU.Triaxial mechanical creep behavior of sandstone[J].Mining Science and Technology,2010,20(3):339－349.

[12]ZHAO YANLIN,CAO PING,WANG WEIJUN,et al.Viscoelasto-plastic rheological experiment under circular increment step load and unload and nonlinear creep model of soft rocks[J].Journal of Central South University of Technology,2009,16(3):484－491.

[13]李亚丽,于怀昌,刘汉东.三轴压缩下粉砂质泥岩蠕变本构模型研究[J].岩土力学,2012,33(7):2035－2040.LI Ya-li,YU Huai-chang,LIU Han-dong.Study of creep constitutive model of silty mudstone under triaxial compression[J].Rock and Soil Mechanics,2012,33(7):2035－2040.

[14]徐进,张家生,赵同顺,等.软弱路基土体三轴蠕变试验及蠕变模型研究[J].中南大学学报(自然科学版),2011,42(10):3136－3142.XU Jin,ZHANG Jia-sheng,ZHAO Tong-shun,et al.Creep triaxial tests and constitutive model of embankment foundation soil[J].Journal of Central South University(Science and Technology),2011,42(10):3136－3142.

[15]熊良宵,杨林德,张尧.岩石的非定常Burgers 模型[J].中南大学学报(自然科学版),2010,41(2):679－684.XIONG Liang-xiao,YANG Lin-de,ZHANG Yao.Non-stationary Burgers model for rock[J].Journal of Central South University(Science and Technology),2010,41(2):679－684.

[16]李娜,曹平,衣永亮,等.分级加卸载下深部岩石流变实验及模型[J].中南大学学报(自然科学版),2011,42(11):3465－3471.LI Na,CAO Ping,YI Yong-liang,et al.Creep properties experiment and model of deep rock with step loading and unloading[J].Journal of Central South University(Science and Technology),2011,42(11):3465－3471.

[17]张志沛,王芝银,彭惠.陕南泥岩三轴压缩蠕变试验及其数值模拟研究[J].水文地质工程地质,2011,38(1):53－58.ZHANG Zhi-pei,WANG Zhi-yin,PENG Hui.Study on triaxial compression creep test and numerical simulation about mud rock in southern Shaanxi[J].Hydrogeology &Engineering Geology,2011,38(1):53－58.

[18]熊良宵,杨林德,张尧.硬岩的复合黏弹塑性流变模型[J].中南大学学报(自然科学版),2010 41(4):1540－1548.XIONG Liang-xiao,YANG Lin-de,ZHANG Yao.Composite viscoelasto-plastic rheological model for hard rock[J].Journal of Central South University(Science and Technology),2010,41(4):1540－1548.

[19]刘姗姗,赵同彬.黏弹性广义Kelvin 模型的FLAC3D二次开发[J].山东科技大学学报(自然科学版),2010,29(4):20－23.LIU Shan-shan,ZHAO Tong-bin.Secondary development on generalized viscoelastic Kelvin model with FLAC3D[J].Journal of Shandong University of Science and Technology(Natural Science),2010,29(4):20－23.

[20]王涛,吕庆,李杨,等.颗粒离散元方法中接触模型的开发[J].岩石力学与工程学报,2009,28(增刊2):4040－4045.WANG Tao,LÜ Qing,LI Yang,et al.Development of contact model in particle discrete element method[J].Chinese Journal of Rock Mechanics and Engineering,2009,28(Supp.2):4040－4045.

[21]张晓平,吴顺川,王思敬.类土质路堑边坡动态监测及数值模拟分析[J].岩石力学与工程学报,2008,27(增刊2):3431－3440.ZHANG Xiao-ping,WU Shun-chuan,WANG Si-jing.Dynamic monitoring and numerical analysis of soil-like cut slope[J].Chinese Journal of Rock Mechanics andEngineering,2008,27(Supp.2):3431－3440.