考虑损伤效应的冻粘土蠕变本构理论研究

文章编号:100121986(2007)0120053203考虑损伤效应的冻粘土蠕变本构理论研究李栋伟,汪仁和(安徽理工大学土木工程系,安徽淮南　232001)摘要:用简单的屈服函数作为冻土塑性损伤变量,对西源模型中的粘塑性蠕变量进行损伤修进。

通过相关推导得到了冻土粘弹塑性损伤本构方程。

通过用户子程序将此本构模型嵌入到大型商业有限元软件ADI NA 中。

有限元计算结果表明:考虑损伤的冻土本构模型能较好的模拟加速蠕变阶段,可用于深井工程实践。

关　键　词:冻土;粘弹塑损伤;本构方程;屈服函数中图分类号:T D265.3 文献标识码:AStudy on frozen soil creep constitutive of theory accounting for damageLI Dong 2wei ,W ANG Ren 2he(Department o f Civil Engineering o f Anhui Univer sity o f Science and Technology ,Huainan 232001,China )Abstract :The simple yield standard is used to determine the damage variables frozen s oil ,therefore X iyuan M odel viscoplastic creep is repaired.The constitutive function of the viscoelastic plastic damage 2coupling is deduced.The viscoelastic plastic damage 2coupling constitutive m odel for geo 2materials is embeded into the commercial FE M s oftware ADI NA as a user ’s subroutine.FE M analysizing suggest that the m odel is of the feature of the non 2linear accelerative creep frozen s oil ,and can be used for engineering practice.K ey w ords :frozen s oil ;viscoelastic plastic damage ;constitutive function ;simple yield function 冻土作为寒区工程与人工冻结工程建设中的一种重要材料,其变形特征的研究已日益受到研究者的高度重视。

自20世纪60年代以来,许多学者致力于这方面的研究,已为工程实践提供了许多可靠的理论依据。

尤其是近年来,随着冻结凿井的深度越来越大,暴露出越来越多的工程问题,牵涉到复杂的地质条件和应力状态,尤其是高应力下的冻土力学问题,原有的概念已不能或完全不能适应工程实际,因此,对冻土力学性能的研究,尤其是冻土本构关系的研究显得非常必要[1]。

吴紫汪[2]等对冻土蠕变变形的细观特征进行了分析;马巍[3]等研究围压对冻土强度特性的影响;中科院冰川冻土研究所何平、程国栋等[4]首次对非饱和冻土的强度特征进行了研究。

目前,国内外虽然对冻土变形特征研究取得了许多成果,但是考虑冻土损伤且用于数值计算的冻土本构还研究较少。

1　冻粘土典型蠕变曲线及蠕变特征本文在自行研制的W3z -200型冻土三轴压缩机上,根据文献[5]进行了相关方法的试验。

图1为冻土典型蠕变曲线:图中蠕变曲线1的应力水平为:σ1=5.5MPa ,σ2=σ3=4MPa ;曲线2的应力水平为:σ1=8.5MPa ,σ2=σ3=5MPa ,由图可以得出如下冻土的蠕变变形规律:图1　冻土典型蠕变曲线Fig.1　T ypical creep curves of frozen s oila .　加载瞬时,都具有一瞬时变形,瞬时弹性应变随应力偏量的增加而增大;b .　当偏应力水平低于某一临界值时,蠕变变形是衰减稳定的,最终蠕变趋于某一稳定值;c .　冻土蠕变在较高的偏应力作用下,会出现加速蠕变阶段,最终式样破坏。

收稿日期:2006207218基金项目:安徽省自然科学基金(050440501)资助作者简介:李栋伟(1978—),男,湖南邵东人,硕士,讲师,主要从事特殊土本构理论和岩土工程数值仿真等教学和研究1第35卷第1期2007年2月 煤田地质与勘探C OA L GE O LOGY &EXP LORATI ON V ol.35N o.1Feb.20072　冻土蠕变损伤变量及考虑损伤的冻土粘弹塑性本构方程2.1　冻土蠕变损伤变量冻土材料服从D -P 屈服准则。

D -P 屈服函数形式如下:F =ασm +3J 2-6c cos φ(3-sin φ),(1)α=6sin φ(3-sin φ),(2)式中　φ为冻土的内摩擦角;c 为冻土的粘聚力,MPa ;σm 和J 2分别为第一应力不变量和第二偏应力不变量。

根据文献[6],利用三轴全应力—应变曲线数据建立冻土粘塑性损伤统计本构关系。

假设冻土微元服从平均概率分布,且初始损伤为0,利用从D -P 屈服准损粘塑性伤变量可以表示为:D P=0 ασm +J 2Φσxασm +J 2-σxσP　σx <ασm +J 2Φσx +σP 1 ασm +J 2>σx +σP,(3)式中　σx 和σP 为平均概率分布的两个参数,根据文献[6]可求得。

2.2　考虑损伤的冻土粘弹塑性本构方程利用流变学中的基本模型,通过组合可以获得很多本构模型来描述岩土的流变特性。

西源模型在岩土工程中应用较多,理论比较完善,如图2。

因此,本文考虑损伤变量时采用西源模型[7-10]为基础,西源模型方程具体形式如下:图2　西源模型Fig.2　X iyuan M odel当σ<σs 时,本构方程为{ε}={σ}E 0+{σ}E 1(1-e E 1t Πη1);(4)当σΕσs 时,本构方程为{ε}=(1E 0+1E 1(1-eE 1t Πη1)+<F >tη2){5Q5σ}。

(5)西源模型从整体上能够描述冻土在高应力下发生的非衰减型蠕变特征,不过由模型理论求得的蠕变本构是线性的,无法描述冻土非线性粘塑性流动性质,也无法描述冻土的细观破坏机理。

因此为了更好的描述特定条件下冻土蠕变特征,需要对蠕变的塑性流动项进行损伤修进。

全解耦合方法认为:损伤对结构中应力、应变场没有影响。

首先利用无损材料的本构关系、平衡方程求解应力场和应变场,将求得的应力和应变代入损伤的演化方程,这种耦合方法最简单。

根据上面推导的损伤变量,并综合西源模型本构方程,采用相关联流动法则,得出冻土的粘弹塑损伤本构方程:{ε}=(1E 0+1E 1(1-e E 1t Πη1){σij }+tη2<F 1-D P >{5F5σij}。

(6)3　用于数值计算的柔度矩阵在实际工程中应用(6)式很难求出理论解。

随着计算机的快速发展,可应用有限元方法来求粘弹塑损伤问题。

而有限元法的核心问题是求解粘弹塑损伤柔度矩阵,因此本文首先介绍粘弹塑损伤柔度矩阵的计算方法。

在计算非线性连续体时,通常把总应变ε分成弹性应变分量εet (包括瞬时弹性分量和粘弹性分量);塑性应变分量εPt (包括瞬时塑性应变分量和粘塑性应变分量)。

因此根据式(6)总应变可表示为:{ε}={εe t }{εPt }。

(7)粘弹柔度矩阵根据文献[11]可以求得:[C e]=1E t 1-μ-μ000-μ1-μ000-μ-μ10000002(1+μ)0000002(1+μ)0000002(1+μ)。

(8)由式(6)得到粘塑性损伤应变的表达式:εPt =<F >t η25Q 5{σ},(9)式中　F 为屈服函数,由式(1)确定,采用D -P 屈服准则,且应用相关联流动法则,粘塑性损伤矩阵采用隐式积分法,可由下面推导求出。

假设在t n时刻,以求得结点位移{δn}和应力{σn},为方便编程采用隐式时间步求解方案,确定Δt =t n +1-t n 时间间隔内产生的应变增量:{Δεn t }P =Δt n {εn +1},(10){εn +1}P =εP t +[C n ]P {Δσn},(11)式中　{Δσn}为Δt 时刻内产生的应力增量,式(11)可写为:{Δεn }=Δt n [{εn }+[C n ]P {Δσn }],(12)式中　[C n]P={5εP5σ}。

(13)・45・煤田地质与勘探 第35卷由式(12)和式(13)可得:[C n ]P=5εP5σ=1η2(5F 5σ)(5F 5σ)T +F 52F 5σ2,(14)式中　5F5{σ}=2sin φ(3-sin φ)[P ]3σm +3[R ]6J 2σ;(15)(5F 5σ)(5F5σ)T =4sin 2φ(3-sin φ)2[P ]9σ2m +14J 2[S ]{σ};(16)52F{σ}2=-2sin φ(3-sinφ)2[P ](9σm )2+3[R ]12J 2J 2{σ};(17)式中　[P ]、[R ]、[S ]为系数矩阵:[P ]=111000111000111000000000000000000000;(18)[R ]=2-1-1000-12-1000-1-12000000600000060000006;(19)[S ]=2-1-1000-12-1000-1-12000000120000001200000012。

(20)由式(3)、(14)得到粘塑性损伤柔度矩阵增量表达式:[C ]P=1η2(1-D P)(5F 5σ)(5F 5σ)T +F 52F 5σ2(21)式中系数矩阵(5F 5σ)(5F 5σ)T 、F 52F5σ2可由式(18)、(19)、(20)求得。

式(8)和式(21)分别表示冻土粘弹性柔度矩阵和粘塑性耦合损伤柔度矩阵。

4　算例及结论对两淮地区深部冻结粘土在三轴恒应力条件下进行了蠕变试验,试验条件轴压:σ1=8MPa ,围压σ2=σ3=5MPa 。

试验结果如图3所示。

冻土本构计算参数如表1所示。

将试验值和理论值比较,如图3所示。

表1　冻土数值计算参数T able 1　Numerical simulation p arameters of frozen soil 试验温度E 0ΠMPaμE 0ΠMPa Η1Π106MPa ・minη2Π108MPa ・minc ΠMPa φΠ(°)-18℃2880.213651.41.73.710图3　冻土三轴蠕变实测值和理论值Fig.3　Observed and theory data of triaxial frozen s oil creep 通过本文分析和数值计算,对复杂应力状态下冻土蠕变损伤本构研究具有如下实用的结论:a .　本文是在冻土宏观蠕变的试验基础上,应用损伤理论与统计学方法,定义了粘弹塑性损伤变量,提出了一个适合冻土材料粘弹塑损伤模型。