基于FLAC3D黄土的湿载结构性本构模型二次开发

罗爱忠， 方 娟

(贵州工程应用技术学院 土木建筑工程学院， 贵州 毕节 551700)

21世纪已经成为土体本构模型特别是结构性土体本构模型研究的世纪[1-5]，针对各种土特性的土体本构模型已经被提出，但是由于没有和计算程序的有机结合，限制了这些模型在实际工程实践中的应用。近年来，为了使现有已开发的本构模型得到更好的应用和被工程界所接受，众多的学者开展了现有本构模型与现有通用软件的有机结合，这其中，又以有限差分法类软件为代表。王春波等[1]推导了硬化土本构模型的有限差分格式，与现有有限差分通用商业软件FLAC相结合进行了二次开发，并将其应用于深基坑稳定性分析。李英杰等[6]考虑了变形模量的劣化，建立了相应的应变软化模型，结合FLAC进行隧道开挖稳定性分析。杨文东等[7]改进了Burgers蠕变损伤模型，并在FLAC软件中实现了二次开发。左双英等[8]在研究现有弹塑性本构模型的基础上，考虑了Zienkiewicz-Pan强度准则，最后基于有限差分法进行了大型地下洞室的开挖稳定性分析。兰航等[9]为了分析节理岩体的采动损伤问题，建立了节理岩体采动损伤本构模型，在FLAC3D中分析了井工开采对上覆岩体的损伤问题。陶慧等[10]将堆石料三维边界面模型嵌入FLAC3D中，进行了堆石料稳定性分析。何利君等[11]将黏弹塑性模型与SMP准则结合，并在FLAC3D中进行了二次开发。陈育民等[12]将邓肯-张模型引入FLAC3D中，并进行了一系列三轴试验验证。综合已有研究，有限差分法由于其概念相对直观，再加上其数学表达式描述得相对简单，在岩土计算领域得到了很大的运用[13-17]。有限差分法直接将所要求解的计算域划分为有限差分网格代替复杂的连续求解域，用网格节点上的函数的差商代替了偏微分方程中复杂的偏导数，从而建立以节点函数值为未知量的代数方程组。湿陷性黄土是一种典型的非饱和结构性黄土，结构性的存在使得其在工程实际中表现出了不同于正常固结土的特殊性质，为了描述黄土结构性的存在对应力应变关系及工程特性的影响，验证证明了该模型在描述黄土应力应变关系方面的准确性和优越性及应用于工程实际分析，有必要将该本构模型与通用计算程序相结合。本文以黄土的湿载结构性本构模型的数值化为切入点，以黄土压剪湿力学作用变化过程分析为目的，验证压黄土的湿载结构性本构模型的合理性和准确性及本构模型与大型数值计算软件FLAC3D结合的可行性与有效性，进而为模型的实际应用奠定基础。

1 FLAC3D分析软件的数值格式及二次开发

1.1 FLAC3D拉格朗日显示差分法

差分法最早由著名学者Wilkins在1963年提出[18]，这一差分格式主要是依据偏导数积分的定义而得到。

在FLAC3D中，单元格式为空间单元，因此其差分格式遵从空间的差分格式，其应变率分量的差分格式可以表述为式(1)形式，FLAC3D中的显式计算中，时间增量步的确定可以通过周期及假设的节点质量确定。

(1)

1.2 FLAC3D分析软件的数值格式

FLAC3D是美国Itasca Consulting Group公司开发的一款用于岩土工程数值模拟的软件。该软件建立在拉格朗日算法的基础上，采用有限差分显式算法来获得模型全部运动方程的时间步长解，该技术和混合离散划分技术的结合，保证了非常精确地模拟塑性破坏和流变，从而追踪材料的渐进破坏，它适合于模拟计算岩土材料力学行为，特别适合模拟材料的高度非线性、不可逆剪切破坏和压密等问题。

FLAC3D中，所有的本构模型都遵从同样的算法。这种算法是在给定t时刻的应力状态和Δt时间步的应变增量，从而求得t+Δt时刻的应力增量和应力状态。假设发生的变形全部是弹性变形，根据胡克定律直接求出应力增量。如果塑性变形发生，根据弹塑性一般原理，应力增量根据总应变增量的弹性部分确定。即，在计算中在每一时步都必须对假设得出的应力增量进行塑性修正，从而得出新的应力增量和应力状态。

(1) 破坏准则。破坏准则满足下式：

f(σn)=0

(2)

式中：f为屈服函数，即当塑性流动发生时的应力组合关系，该函数在应力空间中是一个曲面，所有在曲面之内的应力点保持弹性行为；σn表示σ1，σ2，…，σn。

(2) 应变一般分解为弹性和塑性两部分

(3)

式中：i分别取1、2、3；上标e和p分别表示弹性和塑性。

(3) 应变增量与应力增量之间的弹性关系

(4)

(4) 流动法则。如果采用关联流动法则，及塑性势函数和屈服函数采用相同的形式，则可以得到流动法则满足如下的关系式：

(5)

式中：λ为常数；g为塑性势函数；σ为应力。

(5) 屈服函数的塑性修正。经过塑性修正后，新的屈服函数满足如下的关系：

f(σn+Δσn)=0

(6)

新的塑性应变增量方程被用来评估塑性应变增量的大小。

同时考虑Si的线性关系，得：

(7)

通过流动法则，更进一步应变增量可以表示为：

(8)

式中：Si与应变增量之间的线性关系同样被采用。

硬化常数λ表示为：

(9)

1.3 FLAC3D自定义本构模型的二次开发环境及模型文件编写

在FLAC3D中，可以通过两种方法自定义本构模型，一种方法是采用FISH语言，另外一种方法是采用面向对象的语言C++编写。FISH语言繁琐而且不通用，因而C++就成了FLAC3D中自定义本构模型的主流方法。模型嵌入的方法是将C++编译的模型文件在现有的编译平台上编译成动态链接库文件(DLL文件)，在使用时通过FLAC3D的可执行文件载入。自定义的本构模型的主要功能是给定应变增量，获得新的应力，辅助功能还包括提供模型的名称、版本等基本信息及完成读写的基本操作。

目前FLAC3D中自定义本构模型需要至少在Visual Studio 2005以上的版本进行创建，在FLAC3D中所有的本构模型的加载几乎都是通过动态链接库的形式提供的，采用这种链接方式主要具有以下两个优点：使得用户自定义的本构模型与软件自带的本构模型两者在执行时效率处于同一水平；自定义的本构模型在更高版本的FLAC3D系列软件中得到通用。

2 黄土湿载结构性本构模型在FLAC3D中的计算格式推导

文献[19]提出的考虑压剪湿综合影响的结构性本构模型是基于一定的应力增量得到塑性应变增量，但是在FLAC3D中，其计算格式是基于一定的应变增量得到相应的应力增量。根据弹塑性理论的基本知识，弹性应变增量和塑性应变增量之间的关系可以表示为：

(10)

(11)

依据相关联流动法则，可以得到塑性应变增量表示为：

(12)

由文献[19]的屈服函数，可以得到屈服函数中球应力和偏应力的偏导如下：

(13)

(14)

根据Hooke定律，应力增量与应变增量的关系可以表示为：

(15)

(16)

根据相关联流动法则，新的应力状态可以表示为：

pn=po+Δp

(17)

qn=qo+Δq

(18)

式(18)中，上标n和o分别代表新的应力和塑性修正前的应力。将式(16)和式(17)分别代入式(17)和式(18)得：

pn=pI-Λ*KCa

(19)

qn=qI-Λ*3GCb

(20)

其中，pI和qI分别为：

pI=po+Kdεv

(21)

qI=qo+3Gdεs

(22)

把经过塑性修正后的应力代入屈服面方程，并进一步整理则有：

aΛ2+bΛ+c=0

(23)

其中：

进一步对式(23)求解，可得方程的两个根，若两根都大于零，去最小的一个根，若两根一正一负，取正根。从而，新的应力可以表示为：

(24)

对于湿载条件下的结构性边界面本构模型，由于引入下加载屈服面，不能按常规方程退化求出式(24)，结合相关联流动法则，硬化参量改写为：

(25)

进一步，在主应力空间表示为：

(26)

依据文献[19]湿载条件下的结构性边界面本构模型的屈服函数，有如下的关系式成立：

(27)

式(27)中，各分项满足：

其中，δij为克罗内尔系数。结合式(27)，则有：

式(10)～式(27)中，p，q，σ，ε为土的球应力、偏应力、一般应力和一般应变，下标i和j表示变量顺序；δij为克罗内尔系数；psx为p-q平面上结构屈服面与p轴正半轴交点坐标；pt为p-q平面上结构屈服面与p轴负半轴交点坐标；M为p-q平面上剪切屈服面斜率。

3 黄土的湿载结构性本构模型的实现流程及加载运行

FLAC3D中自定义本构模型的工作空间格式为UDM.dsw，在该工作平台下，以C++编程语言为手段，通过编辑头文件，张量文件及运行文件的塑性指示及流动法则，实现应力应变的计算和应力状态变量的更新。基于上一节的公式推导，结合FLAC3D特点，得到黄土的湿载结构性本构模型的数值实现流程如图1所示。

图1程序实现流程图

模型加载运行过程中，参数的输入必须与程序中一致，见图2。将编好的头文件loess.h和源文件loess.cpp导入到工程文件中，经过编译链接之后，形成了动态链接库文件loesscam.dll，将loesscam.dll拷贝到安装文件目录下，如图3所示。在程序调用时，通过config cppudm进行模块声明，通过model load loesscam.dll语句进行模型加载。判断模型是否被加载成功，可以通过print model 语句判断是否加载成功，如果加载成功将会命令窗口中找到该本构模型及其编号，如图4所示。这样，程序在计算时就可以识别自定义的模型及属性。需要注意的是，在使用RESTORE命令恢复一个含有自定义模型的文件时，也必须使用config cppudm命令和model load loesscam.dll语句。

图2 C++程序中的模型参数

图3 FLAC3D生成动态链接库文件

图4自定义本构模型加载

4 自编本构模型合理性验证

模型的调试与验证实际上包含了两个最主要的方面：一方面是程序编制过程中出现的各种语法及算法的合理性，对于这一方面的问题，参照程序语言的编程规定即可实现；另一方面，需要对自定义本构模型的计算结果进行验证和完善，这个问题的解决方法一般是将模型的计算结果与理论值进行比较，以证明自定义本构模型的合理性和正确性。为了验证自定义本构模型内嵌FLAC3D平台的可行性及合理性，本文采应变式加载三轴试验的方法验证模型开发的正确性。通过FLAC3D外部程序语言Fish语言编写了轴向应力，轴向应变及体应变、剪应变求解的外部函数。模拟的关键在于与黄土常规三轴试验排水条件的相似性。采用的模拟参数见表1，该模型参数通过常规三轴试验求得，鉴于篇幅关系，可参看文献[19]。图5～图8分别给出了含水率10%、固结压力50 kPa和含水率18%、固结压力200 kPa的试验结果及自定义模型模拟结果的对比关系。从图5～图8可以看出，自定义本构模型与试验得到的试验值比较吻合，变化趋势一致，大小相差不大，误差在5%以内，这进一步证明了笔者所建立的黄土的湿载结构性本构模型的合理性及其在FLAC3D平台上实现的可行性和合理性。

表1 计算参数

图5 含水率10%，固结压力50 kPa试验结果

图6 含水率10%，固结压力50 kPa自定义本构模型模拟结果

图7含水率18%，固结压力200 kPa试验结果

参考文献：

[1] 谢定义,姚仰平,党发宁.高等土力学[M].北京：高等教育出版社,2008.

[2] Terzaghi K. Theoretical Soil Mechanics[M]. New York: J. Wiley and Sons, 1943．

图8含水率18%，固结压力200 kPa自定义本构模型模拟结果

[3] Soga K. Mechanical behavior and constitutive modeling of natural structured soils[D]. Berkeley: Univ. of California, Calif. 1994.

[4] Martin D. Liu, Buddhima N. Indraratna, M., General strength criterion for geomaterials including anisotropic effect[J]. International Journal of Geomechanics, 2011,11(3):251-262.

[5] 金福喜，钱 毅，袁权威.某城市排土场边坡三维稳定性分析[J].水资源与水工程学报，2017,28(2):222-226.

[6] 王春波,丁文其,乔亚飞.硬化土本构模型在FLAC～(3D)中的开发及应用[J].岩石力学与工程学报,2014,33(1):199-208.

[7] 李英杰,张顶立,刘保国.考虑变形模量劣化的应变软化模型在FLAC～(3D)中的开发与验证[J].岩土力学,2011,32(S2):647-652，659．

[8] 杨文东,张强勇,张建国,等.基于FLAC～(3D)的改进Burgers蠕变损伤模型的二次开发研究[J].岩土力学,2010,31(6):1956-1964.

[9] 左双英,肖 明,陈俊涛.基于Zienkiewicz-Pande屈服准则的弹塑性本构模型在FLAC～(3D)中的二次开发及应用[J].岩土力学,2011,32(11):3515-3520.

[10] 蓝 航,姚建国,张华兴,等.基于FLAC～(3D)的节理岩体采动损伤本构模型的开发及应用[J].岩石力学与工程学报,2008,27(3):572-579.

[11] 何利军,吴文军,孔令伟.基于FLAC～(3D)含SMP强度准则黏弹塑性模型的二次开发[J].岩土力学,2012,33(5):1549-1556.

[12] 陶 惠,陈育民,肖 杨,等.堆石料三维边界面模型在FLAC～(3D)中的开发与验证[J].岩土力学,2014,35(6):1801-1808.

[13] 陈育民,刘汉龙.邓肯-张本构模型在FLAC～(3D)中的开发与实现[J].岩土力学,2007,28(10):2123-2126.

[14] 邹 博,姚仰平,路德春.变换应力三维化方法在清华模型中的应用[J].岩石力学与工程学报,2005(23):4303-4307.

[15] 杨庆光，李毛毛，王国锋，等．考虑土体结构性损伤的静压桩承载力及时效性研究[J]．湖南工业大学学报，2017,31(6):13-19.

[16] 李忠友，刘元雪，李永毅，等．结构性土在能量转化过程中的损伤演化规律[J]．土木建筑与环境工程，2017,49(5):40-48.

[17] 胡小荣，樊晓梅，童肖龙，等．饱和粘性土的结构性本构模型[J]．南昌大学学报(工学版)，2017,39(3):246-253.

[18] Alder B J, Femback S, Rotenberg M. Fundamental methods in dydrodynamics[J]. Methods in Computational Physics Advances in Research & Applications, 1964:211-263.

[19] 罗爱忠,邵生俊,陈昌禄,等.黄土的湿载结构性本构模型研究[J].岩土力学,2015,36(8):2209-2215.