冻土蠕变模型在ABAQUS中的二次开发

冻土蠕变模型在ABAQUS中的二次开发
曹伟;李文静;张红芬;杜萌洲
【摘要】为了使广义西原模型可以描述冻土的各个变形阶段,用非线性牛顿体替代线性牛顿体进行改进,采用类比的方法将冻土单轴应力状态下的本构方程推广到三维状态;在ABAQUS中利用二次开发平台,编写了改进广义西原模型的UMAT子程序,并在单轴、三轴蠕变条件下进行检验.单轴蠕变的数值解与解析解计算结果十分吻合,两淮地区深部冻结粘土三轴蠕变试验模拟值与实验值相符.表明改进的广义西原模型可以很好地描述冻土蠕变变形特征,包括加速蠕变阶段,UMAT子程序可以用于冻结法施工工程数值模拟.
【期刊名称】《华北科技学院学报》
【年(卷),期】2019(016)001
【总页数】7页(P63-69)
【关键词】广义西原模型;ABAQUS;UMAT;蠕变变形
【作者】曹伟;李文静;张红芬;杜萌洲
【作者单位】华北科技学院建筑工程学院,北京东燕郊065201;燕京理工学院建筑学院,北京东燕郊065201;华北科技学院建筑工程学院,北京东燕郊065201;中国神华能源股份有限公司国华惠州热电分公司广东惠州516001
【正文语种】中文
【中图分类】TU445
0 引言
随着人们对资源的需求急剧增加及冻结法凿井施工技术的日益推广，冻结凿井的深度也越来越大。

冻土中同时有着冰和未冻水，在恒定载荷下冻土发生蠕变，高应力存在会加速蠕变阶段，使冻结管断裂，造成经济损失。

常用的粘弹塑性流变模型有伯格斯模型、Kelvin模型、Bingham模型、西原模型及广义西原模型等。

其中，西原模型和广义西原模型应用较广。

传统的流变模型认为粘滞系数是不变的常数，冻土蠕变的非线性特征、冻土的加速蠕变阶段都不能很好被反映出来，并且与冻土在高应力下实际的蠕变变形相差甚大。

理想的冻土本构模型应考虑时间和应力水平的影响，能较好描述冻土加速蠕变阶段在内的整个流变过程。

因此，对现有流变模型进行改进具有重要的理论意义与工程应用价值[1,2]。

袁文华[3]针对高偏应力下人工冻结粘土的蠕变特性，在西原模型中加入粘弹性元件项，并用抛物线屈服函数替代塑性屈服项进行改进，建立高应力下抛物线型粘弹塑性蠕变本构模型。

汪仁和等[4]结合单轴蠕变试验结果，提出用Mises屈服函数改进西原模型线性牛顿体，描述冻土屈服后的蠕变特性，并依据统计损伤理论建立冻土蠕变损伤本构模型。

姚兆明等[5]将Abel粘壶引入到西原模型，建立分数阶导数西原模型。

李栋伟[6]对标准粘弹塑性模型中的粘塑性项进行改进，建立抛物线型屈服面粘弹塑本构模型，描述冻土在高围压复杂应力状态下的蠕变变形特征，并在ADINA有限元软件中开发了该本构的子程序。

本文对广义西原模型进行改进，采用非线性牛顿体替换线性牛顿体描述冻土的加速蠕变，屈服准则采用Drucker-Prager准则，使之更好地描述冻土的屈服条件。

在有限元软件ABAQUS中，利用二次开发平台，用Fortran语言编写了改进广义西原模型的用户材料子程序(UMAT)，建立了单轴、三轴蠕变压缩试验的模型，并分别与解析解和实验值进行比较，验证了编写的UMAT子程序的正确性。

1 改进广义西原模型理论
1.1 改进广义西原模型
目前，广义西原模型理论较成熟、适用范围较广，但该模型不适用于描述冻土在高应力下的加速阶段的粘塑性变形性质。

为了更好地描述冻土屈服后的蠕变特性，本文采用非线性牛顿体替换线性牛顿体，对广义西原模型进行改进，改进后的广义西原模型如图1所示。

通过大量的冻土试验，可以知道冻土蠕变破坏可以借助Drucker-Prager屈服函数描述[6]。

图1 改进广义西原模型
单轴状态下的改进广义西原模型蠕变本构方程为：
〈〉
(1)
式中EH,E1,E2为弹性模量，MPa；η1,η2为粘弹性粘滞系数,MPa·min；η3也为粘弹性粘滞系数，MPa·min2，〈〉为一个开关函数，当0时，〈〉=0，当>0时，采用类比的方法，一维和三维本构方程间的对应关系为[7-8]: EH↔2GH，
E1↔2G1，E2↔2G2，η1↔2H1，η2↔2H2，η3↔2H3，推得三维应力状态下冻
土的蠕变本构模型为：
(2)
式中GH为瞬时弹性体积剪切模量，MPa；G1，G2为粘弹性体积剪切模量，MPa；H1，H2为粘弹性粘滞系数，MPa·min，H3也为粘弹性粘滞系数，MPa·min2。

1.2 粘弹性和粘塑性应变增量推导
计算非线性连续体变形时通常将总应变增量{Δε}分解成三部分，即瞬时弹性应变增量{Δεe}、粘弹性应变增量{Δεve}和粘塑性增量{Δεvp}之和，如下：
{Δε}={Δεe}+{Δεve}+{Δεvp}
(3)
对于瞬时弹性应变增量有：
{Δεe}=[C]{Δσ}
(4)
式中：[C]为柔度矩阵；{Δσ}为应力增量。

(5)
t+Δt时步Kelvin模型的粘弹性应变{εve}t+Δt可由Kelvin模型的本构关系推得[9]：
(4)
式中：{εve}t为t时步时Kelvin模型的粘弹性应变；Δt为时步增量；{σ}当前时刻的应力[10]。

在Δt内粘弹性应变的增量为:
{Δεve}={Δεve}t+Δt-{Δεve}t
(7)
同样由式(2)得粘塑性应变：
(8)
式中: F为屈服函数，岩土材料F0=1，利用Drucker-Prager屈服准则及相关联流动法则：
(9)
式中:φ为冻土的内摩擦角； c为冻土的粘聚力；σm、J2 分别为第1应力不变量和第2偏应力不变量[11]。

粘塑性蠕变应变增量计算式为：
(10)
由式(9)得：
〈〉
(11)
经推导，
(12)
式中[P]，[R]为系数矩阵：
(13)
(14)
因此，粘塑性应变随时间的增量为
(15)
当粘弹性和粘塑性应变增量确定后，由岩土材料粘弹塑性增量本构关系推得应力增量为：
{Δσ}=[D]({Δε}-{Δεve}-{Δεvp})
(16)
1.3 改进广义西原模型在有限元ABAQUS程序中的实现
有限元软件ABAQUS应用范围广、功能强大，可以模拟由简单线性到复杂非线性的各种问题，当中不仅提供了标准的有限元分析程序，还具有良好的开放性能。

用户可以借助二次开发的平台，开发用户子程序接口和应用程序接口，以此扩展主程序的功能，在实际工程中得到广泛应用[12]。

本文在有限元软件ABAQUS中，利用二次开发平台，用Fortran语言编写了改进广义西原模型的UMAT子程序。

UMAT子程序主要流程如图2所示。

图2 改进广义西原模型UMAT流程图
2 程序正确性验证及算例
2.1 单轴压缩试验算例
为了验证编写改进西原模型UMAT的正确性，将单轴压缩状态下用ABAQUS计算得到的应力应变曲线与解析解进行对比。

建立1 m×1 m×1 m立方体模型，划分一个单元，底端z方向固定，顶面施加分别施加1 MPa、1.5 MPa、2 MPa、2.1 MPa和2.5 MPa的荷载，计算模型的单轴压缩蠕变数值解。

材料参数取值如表1：
表1 材料参数表
GH/MPaνG1/MPaG2/MPaH1/MPa·minH2/MPa·minH3/MPa·min2c/MPaφ/(°)5.00.35.05.00.50.50.51.00
图3 有限元模型及约束与荷载情况
对于蠕变计算，在ABAQUS中设置了两个分析步，类型均为静力分析，第一个分析步总时间设置为1×10-20 min；第二个分析步采用自动增量时间步，初始增量时间步为0.001 min,总时间为1.0 min。

根据材料参数，屈服函数
当材料未屈服，即F<0时，解析解为
当材料屈服，即F≥0时，解析解为
模型的数值解由ABAQUS调用UMAT计算得到，两者的对比如图：
图4 数值解和模型解析解对比
从图4中可以看出，在不同轴压下，数值解和解析解均基本一致。

当轴压小于屈
服应力σs(2 MPa)时，模型处于弹性蠕变阶段，数值解和解析解完全一致，当轴
压大于屈服应力σs(2 MPa)时，模型进入弹塑性蠕变阶段，数值解和解析解略有
差别，但基本一致。

产生差别的原因在于UMAT子程序是基于三维应力状态下的
本构方程建立的，三维应力状态下的本构方程在由一维本构方程类比推导得到时存在一定的近似。

改进的西原模型UMAT子程序很好地描述了材料的瞬时弹性变形，减速蠕变和超过屈服应力后的加速蠕变，数值解与一维的解析解曲线吻合，证明了该子程序在ABAQUS中二次开发的正确性，可以用于工程计算。

2.2 三轴蠕变试验算例
为了验证改进广义西原模型UMAT子程序对于三轴应力状态下的适用性，本文模
拟了三轴蠕变试验，并与对两淮地区深部冻结粘土在三轴条件下蠕变试验结果[13]进行对比。

试验土柱模型为直径36.1 mm，高度80 mm，共划分2000个单元，单元类型选为三维实体单元C3D8，完全积分。

只约束模型底面的z方向的自由度。

三轴蠕变试验的模型如图5所示，数值模拟参数见表2。

图5 有限元模型及网格划分情况
如图6所示，将试验温度-10℃时不同围压和轴压条件下的三轴压缩蠕变试验模拟值与实验值进行对比，可以发现二者在对应围压下曲线吻合度较高。

同时在图可以发现：冻土受偏应力较低时，当前应力状态在屈服面以内，只发生衰减蠕变变形和稳定蠕变变形；偏应力水平较高，当前应力状态超过屈服面时，冻土会产生较大塑
性的流动，出现加速蠕变阶段。

深厚冻土中的地应力往往较大，冻土容易局部进入加速蠕变阶段，故必须考虑。

由三轴蠕变试验模拟值与实验值曲线的对照，进一步论证了采用改进广义西原模型可以很好地描述冻土蠕变变形特征，包括加速蠕变阶段，也验证了UMAT子程序开发的正确性，可以用于工程计算。

3 结论
本文对广义西原模型进行了改进，用非线性牛顿体替代线性牛顿体，并在ABAQUS有限元软件中应用，得到以下结论：
(1) 改进后的模型可以用来描述冻土的每个变形阶段，即衰减蠕变阶段、稳态蠕变阶段和加速蠕变阶段，可以借助此模型模拟分析深井冻结壁的变形规律。

表2 数值模拟参数表
GH/MPaνG1/MPaG2/MPaH1/MPa·minH2/MPa·minH3/MPa·min2c/MPaφ/(°)125.00.28248.0200.04.5e34.5e31.5e62.98.0
(2) 在ABAQUS二次开发平台上，成功地添加了改进广义西原模型的UMAT用户材料子程序。

通过对单轴、三轴蠕变压缩试验的模型的数值模拟，并将数值计算结果分别与解析解和实验值进行了比较，多方面验证了编写的UMAT子程序的正确性，可以应用于冻结法施工工程数值模拟。

(3) UMAT子程序的开发思路可以为其他粘弹塑本构模型的二次开发提供借鉴。

图6 三轴蠕变试验模拟值与实验值对比
参考文献：
【相关文献】
[1] 美启航, 于晖, 陈继.冻土蠕变特性研究现状及展望[J].灾害学, 2018, 33(S1): 47-56.
[2] 马巍.冻结粘性土的变形分析[J].冰川冻土, 2000, 22(1): 43-47.
[3] 袁文华.人工冻土黏弹塑性本构参数反分析研究[J].岩土力学, 2013, 34(11): 3091-3095.
[4] 汪仁和, 李栋伟, 王秀喜.改进的西原模型及其在ADINA程序中的实现[J].岩土力学, 2006,
27(11): 1954-1958.
[5] 姚兆明, 张秋瑾, 牛连僧.基于蚁群算法的冻结重塑黏土分数阶导数西原模型分析[J].长江科学院院报, 2016, 33(07): 81-86.
[6] 李栋伟.高应力下冻土本构关系研究及工程应用[D].安徽: 安徽理工大学, 2005.
[7] 曹伟, 项蒙佳, 赵少飞.冻土蠕变的非线性模型研究[J].华北科技学院学报, 2015,12(3): 82-86.
[8] 郑雨天.岩石力学的弹塑性理论基础[M].北京：煤炭工业出版社, 1988.
[9] 孙钧, 汪炳鉴.地下结构有限元法解析[M].海: 同济大学出版社, 1988: 163-184.
[10] 潘晓明, 杨钊, 雷春娟, 等.广义西原粘弹塑流变模型在ABAQUS中开发与应用[J].建筑结构学报, 2010, 31(S2): 324-329.
[11] 张海银, 范菊红, 李栋伟.动荷载作用下人工冻土蠕变试验研究及有限元分析[J].煤矿安全, 2013, 44(4):195-197.
[12] 曹伟.下加载面修正剑桥模型在ABAQUS中的二次开发及应用[D].哈尔滨: 哈尔滨工业大
学,2014.
[13] 徐大伟.人工冻土黏弹塑性本构模型研究及应用[D].安徽: 安徽理工大学, 2013.。