含分数阶导数元件非线性蠕变模型的二次开发

张 涛，何利军

（南昌航空大学土木与建筑学院，江西 南昌330063）

含分数阶导数元件非线性蠕变模型的二次开发

张 涛，何利军

（南昌航空大学土木与建筑学院，江西 南昌330063）

通过GDS应力三轴仪对湛江黏土在恒定围压不同偏应力下的三轴蠕变试验，得到几组不同偏应力下轴向应变与时间关系数据，鉴于传统线性元件组合模型无法精确描述黏土蠕变的非线性问题，提出采用分数阶导数元件分别替换传统Burger模型中两个牛顿粘壶线性元件，建立起含两个分数阶导数元件的非线性蠕变模型，采用C++语言将该模型在FLAC3D中进行二次开发，通过建立一个三轴数值算例，将不同偏应力下数值计算结果与试验数据进行对比，结果表明提出的非线性蠕变模型相比Burger模型更适合描述湛江黏土的蠕变过程。

分数阶导数；二次开发；蠕变

岩土的蠕变问题一直是工程基础处理、边坡防护和围岩开挖与支护过程中需要重视的棘手难题，而软黏土在工程实践中属于地质不良的土体，其在我国沿海和内陆江河冲积流域分布较为广泛，给工程建设带来诸多潜在危害。

软土的蠕变是岩土流变学研究的重点内容之一，许多研究人员依据大量的试验和工程实践得到了一些可观的成果：如杨爱武[1]对天津滨海吹填软土进行三轴不固结不排水和单轴压缩蠕变试验，结果表明两类试验历时曲线都具有非线性特征；雷华阳[2]对滨海软土进行不同加荷方式下的蠕变试验，得到应变与应力和时间关系，结果表明滨海软土具有明显的非线性特征且其蠕变变形受其结构制约和影响。这些研究成果使得人们对土体蠕变的认识不断加深，在本构模型建立方面也逐渐由单纯的线性理论模型试图描述土体蠕变特性，逐渐发展为半经验半理论模型和非线性理论模型。

传统蠕变模型都是通过粘、弹、塑三类线性元件以串并联方式组合得到，用以描述对应土体的粘弹、粘塑、弹塑等特性，但这些模型具有一定的局限性，只适用于描述线性蠕变过程，而实际软土蠕变具有很强的非线性特性。本文结合前人建立非线性蠕变本构模型研究方法和在FLAC3D中实现二次开发经验，提出采用分数阶导数元件（阿贝尔元件）替换传统Burger模型中的牛顿粘壶，建立起变参数流变本构模型，并在FLAC3D中完成该模型的二次开发工作，通过建立模拟三轴蠕变数值算例，比较不同偏应力下数值结果与试验数据的吻合程度。

湛江黏土是一种具有一定典型性的强结构性黏土[3]，前期通过GDS应力三轴仪对湛江黏土在恒定围压不同偏应力下的三轴蠕变试验，得到几组不同偏应力下轴向应变与时间关系数据[4]，采用各种模型进行了描述[5]，由于传统模型都是线性模型，基于分数阶软体元件构建非线性模型[6]，在理论上具有更好的描述效果。本文在FLAC3D中基于C++进行二次开发，相对于基于FISH语言二次开发而言[7]，其二次开发效果具有更好的稳定性及方便工程技术人员采用的特点。

1 分数阶导数

1.1 分数阶导数简介

分数阶导数理论相比于整数阶导数理论起步较晚。且应用范围并不广泛，在人们长期的科学实践应用过程中，研究人员发现有些问题采用整数阶并不能很好的描述其发展过程，而分数阶能够给出合理的解释且参数简洁明确。直到上世纪80年代，分数阶导数才得到广泛应用，主要集中在天气和气候研究、医学图像处理、工程中非线性问题处理等，很好地解释了自然科学以及工程领域一些非经典现象。由于分数阶导数计算的复杂性及无法向整数阶那样得到精确的解析解，现今对分数阶导数还没有明确的定义，但应用广泛的主要有三种，分别是Riemann-Liouville定义、Grunwald-Letnikov定义和Caputo定义[8-9]，三类定义之间既有联系又有区别。

分数阶相比于传统整数阶理论，优点主要在于：① 能够很好的体现系统函数发展过程；② 与非线性模型相比，物理意义更明确，参数更简洁；③ 克服了经典整数阶理论与试验数据不吻合的严重缺点，使用较少参数即可获得很好的效果[10]。

1.2 分数阶导数元件

分数阶导数元件最主要的特点就是非线性[11-12]，相比传统线性牛顿粘壶元件，更适合描述土体、橡胶、沥青等具有粘性的复杂材料体，该元件最早由挪威数学家阿贝尔提出，又名Abel（阿贝尔）元件，如图1。

图1 阿贝尔元件Fig.1 Abel element

对上述元件，其模型本构关系为

1.3 分数阶导数模型

本文建立的两元件分数阶导数模型，是将Abel元件分别代替经典Burger模型中Maxwell粘壶和Kelvin粘壶而得到（见图2，图3）。

图2 Burger模型Fig.2 Burgermodel

图3 Sburger模型Fig.3 Sburgermodel

对于Burger模型其由Maxwell模型和Kelvin模型串联而成，分别包含一个虎克弹簧体和一个牛顿粘壶，可以用来描述土体的粘弹性，其一维本构方程如下

上式及后式中均统一约定，E表示弹性模量，η表示粘滞系数，G表示剪切模量，下标M和K分别代表Mawell部分参数和Kelvin部分参数

相应三维本构方程如

上式中σ1=σ1-σ3，σ11=σ1+2σ3，其它参数同上。

对于两元件分数阶导数模型本构方程的推导，分两部分，第一部分导出分数阶Maxwell部分本构关系；第二部分导出分数阶Kelvin部分本构关系；上述两部分串联便得到了两元件分数阶导数Burger模型，所以这两部分本构关系的和即为两元件分数阶导数Burger模型本构关系，具体推导如下：

分数阶Maxwell部分由一个虎克弹簧体和一个阿贝尔(Abel)元件串联而成，基于前述Abel元件本构关系，很容易给出该部分一维本构关系如下

式中：ηa为阿贝尔粘壶基础粘滞系数，r为分数阶阶数。

分数阶Kelvin部分由一个虎克弹簧体和一个Abel元件串联而成，其本构关系推导如下[13]

式中：Q＝EK／ηK，D为微分算子，r为分数阶阶数，若r＝1即为整数阶。

将Heaviside单位阶跃函数θ（t）代入式（5）并令θ（t）＝σ（t），表示在t＞0时在模型两端施加恒定的单位应力，通过Laplace变换再反演便得到了式（6）的蠕变柔量表达式

上式中Y1＝1／EK，由展开级数再逐项反演可得，因此蠕变柔量可进一步写为

由于上式含有无穷级数和Gamma函数，形式复杂，不便于直接应用，需对上式进一步简化。在上式无穷级数内添加n=-1项，同时在无穷级数外减去此项，并进一步令r=1-β，0≤β＜1，实际表明软粘土蠕变过程与时间密切相关，为了保留蠕变时间的指数关系，需对无穷级数内的Gamma函数进行简化，令β=0，则上式变为

式（9）即为分数阶Kelvin模型的蠕变柔量，若β＝0即成为线性粘弹性模型，由式（9）和式（4）相加便得到两元件分数阶Burger模型一维本构方程，如下式

相应三维本构方程形式如下

若对式（11）作个简单变换，令ηK＝ηKtβ，ηM=ηaΓ（1＋r）t1-r，便可进一步得到简化后的下式

观察式（12）与Burger三维本构关系式式（3）进行对比发现，经过替换后的Sburger模型本构关系同Burger模型三维本构关系完全一样，两者的区别之一在于分数阶导数模型中粘滞系数ηK＝ηKtβ和ηM=ηaΓ（1＋r）t1-r是随时间成指数形式变化，而Burger模型中这两类粘滞系数始终为常数，基于这一特点，采用C++语言将分数阶导数模型在FLAC3D中实现二次开发成为可能。

2 分数阶导数模型二次开发

2.1 二次开发要点

得出两类模型差异后，将Sburger模型在FLAC3D中进行二次开发，主要工作是通过C++语言建立循环结构[14]，实现参数ηK和ηM随蠕变时间呈指数关系变化。借鉴FLAC3D中内置的Burger模型源程序(包含. cpp源文件和.h头文件)编写规则和各函数功能，以该源程序为蓝本通过if..else条件选择结构和for循环语句实现两类粘滞系数随时间的变化关系。并将更改后的分数阶模型程序文件以项目形式添加入visual studio 2008中生成Dll(动态链接库)文件[15]，并命名为userSburger.dll。然后将该链接复制到FLAC3D内置模型程序文件目录下，通过在FLAC3D窗口输入相关命令，完成所开发模型同FLAC3D的接口、载入，初步实现该模型二次开发。

模型开发过程中对程序的更改主要集中在基本参数的定义、常量与过程量的定义、基类的描述、模型注册ID号和Run（）函数中本构模型计算关系的修改等。

2.2 二次开发验证

在开发模型同FLAC3D进行接口、载入后，需要对模型开发的正确性及用于数值计算的可靠性进行验证。验证方法分为两部分，一是判断所开发模型程序编写是否正确，即该模型能否用于数值计算；其次判断所开发模型是否实现两类粘滞系数随时间成指数形式变化，这是开发成功的关键。首先通过建立一个简单算例，在Sburger模型中对两个分数阶阶数γ和β取值0退化成Burger模型形式，得ηK=ηK和ηM=ηa；保持其它参数不变，计算结果如图4。

图4 Sburger模型退化后计算结果Fig.4 Calculation of deteriorated Sburgermodel

从图4可以看出，对Sburger模型的分数阶阶数取值为0后，保持其他参数相同，其数值计算结果与内置Burger模型完全吻合。说明两类分数阶模型按照上述方式通过取特殊值可退化为经典Burger模型，也进一步表明开发的两类分数阶模型是可以在FLAC3D中作为本构模型来选取并参与数值计算。

图5 Sburger模型中不同阶数数值计算曲线Fig.5 Numerical calculation curves of different orders in Sburgermodel

从图5可得：对Sburger模型阶数取不同值，数值计算曲线存在明显差异，其基本规律为当β取小，γ取大时，试样在等速蠕变阶段的变形速率较大；而β取大，γ取小时，试样在等速蠕变阶段的变形速率较小；出现这一规律的原因在于β取值越小使得ηK随时间增加的变化率下降，即单位时间内的粘滞系数增长率降低，导致在恒定偏应力作用下，其对应蠕变变形量增加；反之β取值越大，单位时间内粘滞系数增长率较高，使得对应蠕变变形量降低。同理对于另一阶数γ，其取相同值所得ηM的结果正好与β相反。综上，两个分数阶阶数的不同取值，将影响两类粘滞系数随时间变化增长速率的快慢，从而在恒定偏应力下，使得等速蠕变阶段蠕变变形存在显著差异。

3 三轴数值算例

3.1 三轴蠕变试验数据

两类模型的基本参数均依据三轴蠕变试验数据通过各模型三维本构关系拟合得到，表1为湛江软粘土试样基本物理力学性质[16]，图6为该试样在50 kPa恒定围压不同偏应力下，柱体轴向应变量随时间变化关系。

表1 湛江软粘土试样基本物理力学性质Tab.1 Basic physical andmechanical properties of Zhanjiang clay samp le

图6 三轴蠕变试验曲线Fig.6 Curves of triaxial creep test

图6是对试样进行分级加载条件下得到，而只有分别加载才能得到更加真实的土体蠕变曲线，且考虑到岩土材料流变的非线性特性。需将分级加载下的蠕变数据转化为分别加载下的蠕变数据，这里采用陈宗基先生提出并由其学生发展的“陈氏法”[17]进行处理，整理后的蠕变数据如图7。

图7 经“陈氏法”整理后三轴蠕变试验曲线Fig.7 Curves of triaxial creep test under Chen’smethod

3.2 模型建立

为了对比分数阶导数模型与Burger模型描述湛江软粘土蠕变过程的可靠性，在FLAC3D中建立一个圆柱体模型(如图8)算例模拟三轴蠕变过程。圆柱高为2个单位，长和宽均为0．5个单位，半径为0．25个单位，先对模型施加50 kPa的围压和重力加速度，然后模型各方向位移清零，再在上顶面施加竖直方向压力，蠕变时步取为1e-4，分别采用Burger模型和本文开发的分数阶导数Sburger模型进行计算。模型参数取值依据相应本构关系基于三轴蠕变试验数据拟合得到。

图8 柱体模型示意图Fig.8 Schematic diagram of cylindermodel

3.3 参数选取

采用1Stopt数学软件及孙均[18]提出的拟合方法拟合得到两类模型的基本参数，大致可分为两部分，第一部分当t→o时，只有Maxwell弹簧存在变形，可拟合出参数K和GM，第二部分当t→∞时，蠕变达到稳定，将此时本构关系式减去原本构关系式便得到包含粘性和粘弹性本构关系式，与对应的三轴试验数据进行回归拟合。得到Burger模型在不同偏应力下参数如表2，同理可得分数阶Sburger模型参数如表3。并由拟合相关系数R可以看出，各参数拟合情况较好。

表2 恒定围压不同偏应力下burger模型相关参数Tab.2 Parameters of Burgermodel under different partial stress of constant confining pressure

表3 恒定围压不同偏应力下Sburger模型相关参数Tab.3 Parameters of Sburgermodel under different partial stress of constant confining pressure

由表2和表3中两个模型的相关拟合系数R可知，Sburger模型的拟合效果要优于Burger模型。其原因在于湛江黏土的蠕变包含有非线性变形部分，而传统线性模型无法准确描述非线性问题，所以采用非线性本构关系拟合参数的效果要比线性本构关系好；其次，Burger模型在30 kPa偏应力下，拟合系数较高，这是在低偏应力下，黏土蠕变主要产生线性弹性、粘性和粘弹性变形；而当偏应力增加，随着蠕变变形的进一步增长，非线性变形也在不断上升成为后期变形的主要增长部分，是导致65～135 kPa偏应力下Burger模型拟合效果较低的原因之一。另外，从参数拟合效果来看，在数学角度上能一定程度说明Sburger模型本构关系要比经典Burger模型更能描述试验曲线。

将上述两模型采用对应参数分别在不同偏应力下进行三轴数值计算，并将两者数值计算结果同三轴蠕变试验数据进行对比，如图9～图12。

由图9～图12不难看出，开发两元件分数阶导数模型比传统经典Burger模型，更适合描述湛江软粘土的初始蠕变和稳定蠕变阶段。从而验证了本文提出采用分数阶元件替换经典整数阶模型中牛顿粘壶元件，建立变参数流变模型，描述软土非线性蠕变思路的正确性；也进一步说明，本文对分数阶导数模型在

FLAC3D中实现二次开发的正确性和可行性。

图9 偏应力30 kPa下数值计算曲线Fig.9 Numerical calculation curve under the partial stress 30 kPa

图10 偏应力65 kPa下数值计算曲线Fig.10 Numerical calculation curve under the partial stress 65 kPa

图11 偏应力100 kPa下数值计算曲线Fig.11 Numerical calculation curve under the partial stress 100 kPa

图12 偏应力135 kPa下数值计算曲线Fig.12 Numerical calculation curve under the partial stress 135 kPa

4 结论

从本文提出分数阶元件替换传统线性粘壶元件，建立变参数蠕变模型，并在FLAC3D中完成二次开发和验证的结果来看，分数阶元件相对于传统元件，对蠕变曲线的刻画能力更强，预示着对工程的数值模拟结果更加准确。由于土体蠕变变形具有广泛的非线性特性，分数阶元件相对于整数阶元件，将更具备普遍意义。伴随着工程建设的需要和蠕变模型自身的发展，本文研究成果将在未来的岩土工程设计计算中被越来越多地接纳和采用。

[1]杨爱武，张兆杰，孔令伟.不同应力状态下软黏土蠕变特性试验研究[J].岩土力学，2014，35（S2）：53-60.

[2]雷华阳，贾亚芳，李肖.滨海软土非线性蠕变特性的试验研究[J].岩石力学与工程学报，2013，31（1）：2806-2816.

[3]拓勇飞，孔令伟，郭爱国，等.湛江地区结构性软土的赋存规律及其工程特性[J].岩土力学，2004（12）：1879-1884.

[4]蔡羽，孔令伟，郭爱国，等.剪应变率对湛江强结构性黏土力学性状的影响[J].岩土力学，2006（8）：1235-1240.

[5]何利军，蔡羽，吴应红.湛江软黏土的流变本构模型研究[J].合肥工业大学学报：自然科学版，2011（4）：565-570.

[6]何利军，孔令伟，吴文军，张先伟，蔡羽.采用分数阶导数描述软黏土蠕变的模型[J].岩土力学，2011（S2）：239-243.

[7]何利军.分数阶导数蠕变模型及在FLAC3D中的数值实现[J].南昌航空大学学报：自然科学版，2015（4）：35-39.

[8]李秀.几类粘弹性材料的分数阶积分模型及其数值解法[D].宁夏：宁夏大学，2013：8-24.

[9]王玉娇.分数阶导数及其应用[D].太原：太原理工大学，2014：10-44.

[10]陈文，孙洪广，李西成.力学与工程问题的分数阶导数建模[M].北京：科学出版社，2010：124-140.

[11]郭佳奇，乔春生，徐冲，等.基于分数阶微积分的Kelvin-Voigt流变模型[J].中国铁道科学，2009，30（4）：1-6.

[12]何志磊，朱珍德，朱明礼，等.基于分数阶导数的非定常蠕变本构模型研究[J].岩土力学，2016，37（3）：737-744.

[13]王志方，张国忠，刘刚.胶凝原油的黏弹性流变模型[J].高校化学工程学报，2008，22（2）：351-355.

[14]古万荣著.Visual C++完全自学手册[M].北京：机械工业出版社，2009：29-111.

[15]陈育民，徐鼎平.FLAC/FLAC3D基础与工程实例[M].2版.北京：中国水利水电出版社，2013：50-165.

[16]孔令伟，何利军，张先伟.湛江黏土的蠕变模型与变参数塑性元件[J].岩土力学，2012，33（8）：2241-2246.

[17]陈宗基，康文法.岩石的封闭应力、蠕变和扩容及本构方程[J].岩石力学与工程学报，1991，10（4）：299-212.

[18]孙钧著.岩土材料流变及其工程应用[M].北京：中国建筑工业出版社，1999：146-158.

Secondary Development of Non-linear Creep Modelw ith Fractional Order Derivative Elements

Zhang Tao，He Lijun
（School of Civil Engineering and Architecture，Nanchang Hangkong University，Nanchang 330063，China）

Triaxial creep test of Zhanjiang clay was conducted at constant confining pressure and under different deviatoric stress through GDS,which obtained several groups of data concerning the relationship between axial strain and time under different deviatoric stress.In view of the fact that traditional linear elements of nonlinear creepmodel does not accurately describe the clay problem,a fractional order derivative elementwas proposed to replace the two linear elements of the Newton’s pot in the Burgermodel respectively.Then a non-linear creepmodel with two fractional derivatives was established.By using C++language,the secondary development of themodel in the FLAC3D was completed.Finally,comparison wasmade between numerical results under different stress and experimental data by establishing a triaxial numerical example.The results show that compared with the Burgermodel,the non-linear creepmodel ismore favorable for the description of creep process for Zhanjiang clay.

fractional derivative；secondary development；creep

1005-0523（2017）05-0021-08

T U433

A

2017－05－09

南昌航空大学校级基金项目（EA201111317）

张涛（1991一），男，硕士研究生，主要研究方向为蠕变模型的二次开发。

何利军（1977一），男，讲师，博士，研究方向为结构性土体本构模型研究。

（责任编辑 王建华）